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ABSTRACT 


This  thesis  reviews  various  algorithms  for  computing 
solutions  of  systems  of  linear  algebraic  equations  and 
for  computing  the  inverse  of  general  matrices,  block- 
symmetric  matrices,  and  band  matrices.  Norm  bounds  for 
the  round-off  error  incurred  in  performing  these  computa¬ 
tions  are  obtained.  Inverses  of  several  matrices  and 
solutions  of  several  systems  of  equations,  are  computed 
numerically,  and  a  comparison  is  made  between  the  predicted 
bounds  of  error  and  the  actual  round-off  errors. 
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CHAPTER  I 


INTRODUCTION 


1 . 1  History 

The  development  of  automatic  digital  computers  has 
tremendously  increased  the  capability  of  a  numerical 
analyst  to  carry  out  numerical  solution  of  a  mathemati¬ 
cal  problem  involving  a  very  large  number  of  arithmetic 
operations.  But,  due  to  the  amount  of  computation 
involved,  the  intermediate  results  are  so  numerous  so 
as  to  preclude  careful  scrutiny  by  the  analyst.  Conse¬ 
quently,  it  could  happen  that  due  to  the  accumulation 
of  round-off  errors,  the  final  results  obtained  may  be 
completely  meaningless  in  terms  of  the  original  problem. 
Therefore,  if  the  speed  of  an  electronic  digital  computer 
has  relegated  the  simplicity  of  the  numerical  computation 
to  a  somewhat  secondary  position,  it  has  also  established 
a  need  for  studying  the  cumulative  effect  of  round-off 
errors  made  in  performing  numerical  computations. 

The  problem  of  estimation  of  round-off  errors  in 
the  solution  of  linear  algebraic  equations  and  inversion 
of  matrices  was  first  considered  by  Von  Neumann  and 
Goldstine  [13],  and  Turing  [12].  The  former  gave  rigorous 


- 


> 


- 


' 


* 


■ 


. 

. 


s. 

. 

. 

. 


2 


error  analyses  of  the  techniques  based  on  Gaussian  elimina¬ 
tion  and  may  be  said  to  have  laid  the  foundation  of  modern 
error  analysis.  An  outstanding  contribution  to  the  field 
of  error  analysis  has  been  made  by  Wilkinson  [14-19]  whose 
researches  on  this  topic  have  been  the  subject  matter  of 
several  papers  and  books . 

1  c  2  Purpose  of  Study 

An  important  criterion  for  the  selection  of  a  parti¬ 
cular  method  for  the  solution  of  a  system  of  linear 
algebraic  equations  and  inversion  of  matrices  is  that  the 
method  which  leads  to  smaller  error  bounds  for  the  round¬ 
off  errors  incurred  in  performing  numerical  computations 
is  better,  ot-her  considerations  being  equal.  However,  the 
actual  errors  depend  on  the  nature  of  elements  of  the 
matrices  considered.  In  this  paper,  we  attempt  to  find 
upper  bounds  for  the  round-off  errors  incurred  in  the 
computation  of  the  solution  of  a  system  of  linear  algebraic 
equations  and  inversion  of  (a)  general  matrices  using 
Schur's  identity  [ 9  J ;  (b)  block-symmetric  matrices  using 

computational  procedures  due  to  Charmonman  [1];  (c)  band 

matrices  using  minimum  storage. 
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CHAPTER  II 


BASIC  THEORY 

2 . 1  Matrix  Algebra 

We  refer  in  this  section,  to  some  conventions  of 
notations  and  discuss  some  of  the  properties  of  matrices 
which  will  be  required  repeatedly  in  the  remainder  of 
this  paper.  Some  important  definitions  including  the 
concept  of  matrix  norm  of  square  and  rectangular  matrices 
are  also  introduced. 


2.1,1  Basic  Notations  and  Definitions  We  denote 

matrices  by  upper-case  Roman  and  Greek  letters,  and 

vectors  and  scalars  by  lower-case  Roman  letters.  The 

(i,j)  element  of  a  matrix  A  is  denoted  by  a...  We 

—  J 

say  that  a  rectangular  matrix  A  =  (a^,)  consisting  of 
m  rows  and  n  columns  is  of  dimension  m-by-n.  If  the 
number  of  rows  and  columns  of  the  matrix  are  equal,  the 
matrix  A  is  square  and  of  order  n.  In  general,  matrices 
will  be  square  and  of  order  n,  and  vectors  of  dimension 
n,  unless  it  is  indicated  otherwise.  The  identity  matrix 
is  denoted  by  I  and  the  null  matrix  by  0.  The  nota¬ 


tion 


A 


is  used  for  the  matrix  which  has 


a.  .  . 

-ini’ 


the 


absolute  value  of  a. .  for  its  (i,j)  element  and  the 

ij  30 


' 


; 


' 

; 


■ 

■ 
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notation  det  (A)  for  the  determinant  of  a  square 

matrix  A.  The  inverse  of  a  non-singular  square  matrix 

A  is  denoted  by  A  ^ .  If  the  matrix  A  is  real,  then 
T 

A  denotes  the  transpose  of  A,  The  matrix  A  is  called 

T 

symmetric  when  A  =  A  .  A  matrix  A  of  order  2n  is 

called  block-symmetric  when  it  can  be  partitioned  as 
C  !  D 

g-j”  ,  where  C  and  D  are  square  matrices  of  order 

n.  A  matrix  A  is  said  to  be  a  band  matrix  of  band-width 

(2p+l)  if  a.  o  =0  for  all  i  and  j  such  that 

J 

li-J^P"  We  denote  by  L  a  lower-triangular  matrix  in 
which  all  elements  above  the  main  diagonal  are  zero.  The 
lower-triangular  matrix  L  is  called  unit  lower-triangular 
if  the  elements  on  the  main  diagonal  are  all  unity. 
Similarly,  we  denote  by  U  an  upper-triangular  matrix  in 
which  all  elements  below  the  main  diagonal  are  zero,  and 
if  the  main  diagonal  elements  are  all  unity,  the  matrix 
U  is  called  unit  upper-triangular. 

2.1.2  Norms  of  Vectors  A  norm  of  a  vector  x  is  a 
non-negative  scalar  ||xj|  which  is  a  measure  of  the  size 
of  the  vector,  and  satisfies  the  conditions 


. 
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(i)  ||  x  |  >  0  for  x  ^  0  and  jjo||  =  0; 

(2.1)  (il)  ||  k  x  j  =  |k|  j|  x  ||  for  any  scalar  k; 

(ill)  ||  x+y  J  <  |  x||  +  ||  y  |  . 

Three  different  norms  of  vectors  which  satisfy  the 
above  conditions  are : 

(2-2)  1.x |  =  l  |x  |  3 

1  =  1 

(2.3)  II  x  ||  ^  =  max  |x±  |  , 


and 

(2.4)  NIe  =  (  l  |x1|2)^  . 

1=1 

2.1,3  Norms  of  Matrices  A  norm  of  a  square  matrix 
A  is  a  non-negative  scalar  j|  A ||  which  is  a  measure  of 
the  magnitude  of  the  matrix 3  and  satisfies  the  conditions 


(i)  || A ||  >  0  if 


A  ¥■  0  and  ||0|j  =  0  ; 


■ 

■ 

■ 

' 


1 1  -  l 


. 


■ 


(2.5)  (ii)  II  kA||  =  | k  |  ||  A||  for  any  scalar  k  ; 

(ill)  ||A+B||  <  II  A|  +||  B||  ; 

(iv)  ||  AB ||  <  ||  A||  ||  B ||  . 


The  following  matrix  norms  are  in  common  use, 
each  of  which  satisfies  the  above  conditions: 


(2.6) 

II  A||  , 

n 

=  max  l 

J  i=l 

|  a .  .  1  , 

1  ij  1  ’ 

(2.7) 

II  AIL 

n 

=  max  l 

i  ,1  =  1 

| a.  .  |  , 

1  ij  1  5 

(2.8) 

f\l 

=  (maximum 

T  t' 

eigenvalue  of  A  A) 2  , 

and 

(2.9)  II  All  =  (  l  l  |a  |2h  . 

^  1=1  j=l  1J 

The  definition  of  matrix  norm  has  been  extended 
by  Householder  [6],  and  Korganoff  and  Pavel-Parvu  [8]  t 
include  rectangular  matrices .  A  rectangular  matrix 


can  be  normed  by  adjoining  null  rows  and  columns  so  as 
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to  make  it  a  square  matrix.  Once  this  is  done,  the 
definition  of  norm  of  a  square  matrix  may  be  applied, 

2.1,4  Some  Elementary  Properties  of  Norms  In  this 
section  we  assemble  certain  inequalities  involving  vector 
and  matrix  norms  which  we  will  use  repeatedly.  Through¬ 
out  this  paper,  we  will  only  use  the  row  norm  or  “-norm, 
since  it  can  be  easily  computed.  The  following  properties 
hold  for  norms  defined  in  Subsection  2,1,3,  unless  it  is 
indicated  otherwise.  For  the  proof  of  these  properties 
see  for  example,  [11]  or  [16], 


(A) 


(B) 


(C) 


and 


||  |A|j|  =  ||  A||  ,  for  any  of  our  norms 


except  the  spectral  norm,  (2.8), 


■"  '  I  ’  \ 


' 


' 


' 


- 
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(D) 


If  ||Aj|  <  1,  then  (I+A)  is  non-singular. 


(E) 


If  ||  A ||  <  1,  then 


(i) 


||  (I+A) 


-i"  <  _L 


l- 


and 


(ii) 


1  ( I+A)  1  -  l||  < 


Ml 


1- 


(F)  If  A  is  non-singular  and  j|A  ^||  ||E||  <  1  , 

then  ( A+E )  is  non-singular  and  hence 


II  A  1 II II  A  1j|  ||  E|| 

i-I|a“1||||e|| 


(i)  II  (A+E)  1  -  A  1 1|  < 


and 


(ii) 


(A+E) 


II  a"1!! 


< 


i— 1|  a  1||!|e!| 


2 . 2  Gaussian  Elimination 

In  this  section  we  describe  one  of  the  best  known 
and  most  widely  used  'direct  method’  for  solving  linear 
systems  of  algebraic  equations  and  for  inverting  matrices. 
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due  to  Gauss  [5].  A  direct  method  is  one  in  which  after 
a  finite  number  of  operations  a  computed  result  is 
obtained  which  would  be  the  same  as  the  true  solution 
of  the  given  system  if  the  computations  are  carried 
without  round-off.  Basically,  the  Gaussian  elimination 
method  is  simple.  The  first  equation  in  the  given 
system  is  used  to  eliminate  the  first  unknown  from  the 
last  (n-1)  equations,  then  the  new  second  equation  is 
used  to  eliminate  the  second  unknown  from  the  last  (n-2) 
equations,  etc.  The  process  of  forward  elimination  is 
completed  in  (n-1)  reductions.  The  resulting  system 
which  is  equivalent  to  the  original  system  is  upper- 
triangular  and  can  be  solved  by  back-substitution. 


Let 


a(i) 

all 

a(1; 

a12 

a(1) 

o  0  •  CX  -s 

In 

(2.10) 

11 

1 — 1 

< 

11 

< 

aU) 

a21 

0 

a':i) 

a22 

0 

a(1) 

poo  CX 

2n 

0 

0 

(1) 
a  , 
nl 

_ 

a 

a<£> 

e 

0 

(i) 

.  .  .  a 

nn 

be  the 

coefficient  matrix 

of  a 

system  of 

of  order  n  and  rank  n.  Then  after  k 


* 


■ 

-J 


linear  equations 
reductions ,  the 


■ 

: 

. 
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reduced  matrix  A^k^ 


(2.11)  A(k)= 

0 

« 

0 


0 


0 


becomes 

a(1) 

a12 

a(1) 

alk-l 

a(1) 

lk 

a(2) 

22  '  0  ' 

a(2) 

2k-l 

a(2) 

a2k 

a(2) 

2n 

0 

c 

• 

O 

0 

0 

A 

• 

♦ 

0 

(k-1) 
k-1, k-1 

A 

a(k-l) 

ak-l,k' ‘ ’ 

• 

a(k" 

k-1 

0 

• 

0 

• 

• 

• 

• 

• 

0 

0  cep 

e 

0 

a(k> 

nk 

a<k> 

nn 

where 


(2,12) 


(k-1) 

ij 


(k-1) 

ik 


(k-1) 

kk 


(k-1) 

ij 


k=  2  3 


j=k} 

l=k. 


3n  3 

,n  . 


The  solution  x  =  ( x^ , x^ , . » , x  )  to  the  system  of  equations 
Ax  =  b  is  easily  computed  from  the  n1^  reduced  matrix 
A(n)  by 


(2.13) 


x . 

l 


1 


11 


[a 


(1) 

1  j  n+1 


n 

I 

j=i  +  l 


al1 jx.] 

iJ  J J 


i=n, 


...  j  1  3 


where 


ai^n+la  j • «  •  ,n  denote  the  components  of 


- 

. 

i 

. 


- 
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the  right-hand  vector  in  the  n^  reduced  system. 

The  algebraic  basis  of  the  Gaussian  elimination  is 
contained  in  the  following  theorem: 


Theorem  2 . 1  Given  a  square  matrix  A  of  order 
ftj  let  A^  denote  the  principal  minor  matrix  constructed 
from  first  k  rows  and  columns .  Assume  that 


det(A^)^0  for  k=  1  ,  2  ,  «,  .  .  ,  n- 1 

unit  lower- triangular  matrix 

triangular  matrix  U=(u..)  s 
n 

det(A)  =  tt  w...  where  det 
.  -  vv3 

V=1 

of  matrix  A. 

A  proof  of  this  theorem 


.  Then  there  exists  a  unique 
L=(m.o)  and  a  unique  upper- 
o  that  LU=A.  Moreover 3 
(A)  denotes  the  determinant 

is  given  in  [4 ] . 


2.2.1  A  Need  for  Pivoting  In  using  Gaussian 
elimination  with  natural  ordering  we  have  Ignored  the 
possibility  of  the  break-down  of  the  process  due  to  one 
of  the  diagonal  elements  having  the  value  zero.  Ignoring 
the  possibility  of  vanishing  of  the  last  diagonal  element 
in  the  nth  reduced  matrix  A''n\  in  which  case  the  matrix 
is  singular,  let  us  consider  that  an  intermediate  diagonal 
term  becomes  zero.  It  means  only  that  a  leading  submatrix 


. 

. 


r.  '■ 


' 


, 

■ 


-  : -■  y :  ■ 
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of  A  has  a  zero  determinant,  and  does  not  necessarily 
imply  the  singularity  of  A.  Let  us  assume  that 
aH=0  and  since  det(A)^0,  we  can  find  an  a^-^0  for 
some  i>l.  Thus  zero  pivotal  element  can  be  avoided  by 
interchange  of  the  i ^  row  with  the  first  row.  This 
procedure  can  be  employed  at  all  stages  of  reduction. 

In  theory,  interchange  is  necessary  only  when  a  pivotal 
element  is  exactly  zero.  However,  In  practice,  this 
procedure  is  satisfactory  only  when  the  magnitude  of  the 
elements  of  A  and  the  successive  derived  systems  do 
not  vary  appreciably  in  size.  But,  when  any  of  the 
pivotal  elements  ,  afjP  is  close  to  zero  or  small  in 
magnitude  compared  to  a^. ^  ,  k>i,  then  a  serious  round¬ 
off  error  may  be  incurred  since  the  percentage  error 
in  the  reciprocal  of  a  number  varies  inversely  with  the 
magnitude  of  the  number.  To  avoid  this  we  may  use 
pivoting  strategies,  either  partial  or  complete,  [18], 

2.2. 1.1  Partial  Pivoting  At  the  ith  stage 
in  the  reduction  of  A  Into  a''^,  we  locate  the  element 
of  maximum  magnitude  in  rows  i  through  n.  Suppose  the 
pivotal  element  is  in  the  jth  row,  then  rows  i  and  j 
are  inter changed .  The  interchange  of  rows  can  be 


- 

* 


■ 


. 


( 
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accomplished  by  premultiplying  the  i^  reduced  matrix 

by  the  permutation  matrix  FL  ^ .  However,  in  practice 

physical  interchange  of  rows  can  be  avoided  by  constructing 

a  vector  of  integers  p^ ( i=l , 2 , . . ,n)  and  calling  for 

a^  .  instead  of  a. . . 

PiaJ  ij 

2. 2. 1.2  Complete  Pivoting  By  complete  pivoting 

t  h 

we  mean  the  following:  at  the  i  stage  in  the  reduction 

( k ) 

of  A  into  A  ,  we  locate  the  element  of  the  maximum 

magnitude  in  rows  i  through  n  and  columns  i  through 

n,  i.e.,  we  find  the  element  of  maximum  modulus  in  the 

matrix  formed  by  .deleting  the  first  i-1  rows  and  columns 
/  •  \ 

of  A'  .  The  pivotal  element  is  transferred  to  position 

/  •  \ 

(i,i)  of  matrix  A^  '  by  one  row  interchange  and  one 
column  interchange.  This  interchange  can  be  effected  again 
by  use  of  permutation  matrices. 

2 . 3  Error  Analysis 

In  this  section  we  give  an  account  of  the  error 
analysis  of  some  of  the  fundamental  mathematical  computa¬ 
tions  performed  on  a  digital  computer  using  floating¬ 
point  arithmetic.  We  consider  only  floating-point  arith¬ 
metic  since  it  is  available  in  all  modern  digital 


'  ■  ■ 

. 


V 


. 
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■ 
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computers.  Moreover,  fixed-point  arithmetic  is  not  suita¬ 
ble  for  handling  numbers  of  large  and  varying  magnitudes 
in  scientific  computation.  All  the  results  presented 
in  this  section  are  due  to  Wilkinson  [17]. 

There  are  two  main  forms  of  error  analyses  that  can 
be  used.  We  may  attempt  to  trace  the  forward  propogation 
of  individual  rounding  errors  and  then  compare  the  computed 
numbers  with  those  which  exact  computation  would  produce 
or  we  may  show  that  computed  results  may  be  obtained  by 
performing  the  exact  computation  on  a  perturbed  problem. 

The  former  is  referred  to  as  ’forward  analysis'  and 
the  latter  'backward  analysis'.  We  will  only  consider 
the  backward  error  analysis  since  this  is  often  much 
simpler,  particularly  if  floating-point  arithmetic  is 
used.  Assume  that  we  calculate  a  quantity  x  given  by 
the  mathematical  expression 

(2.14)  x  =  f(Xl,x2,...,xn)  . 

It  follows  using  backward  analysis  that  the  computed  x 
satisfies  exactly  an  equation  of  the  form 

=  f  (Xj+e^Xj+Ej  , .  .  .  >xn+en) 


(2.15) 


X 


3 


■ 


' 

'  ■ 


' 


' 


’  :•*  ,•  ff  L. ‘lr‘  '  >'  fc  r.  -19':  on  3.L  ^  -  -V 


' 

’ 
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where  are  the  perturbations.  We  attempt  to  place 

bounds  on  the  magnitude  of  the  perturbations  • 

Corresponding  to  the  mathematical  equation  (2.14) 
defining  x,  we  have  a  computational  equation  which 
defines  the  computed  value  of  x,  as  in  (2.15)  followed 
by  bounds  for  ei ( i=l , 2 , . . ,n) . 

In  binary  floating-point  arithmetic  each  number  x 
is  represented  in  the  form  x=2°(a),  where  b  is  a 
positive  or  negative  integer,  and  a  is  a  number  which 
satisfies  h  <  |a|  <  1.  The  number  a  is  called  the 
mantissa  or  fractional  part  and  b  is  called  the 
exponent  or  index,  of  x.  We  will  denote  by  t,  the 
number  of  binary  digits  in  the  mantissa  of  a  floating¬ 
point  number. 

2.3.1  Common  Floating-Point  Operations  Let  the 
number  e  denote  the  unit  round-off  error,  i.e., 

(2.16)  | e |  <  2-t  , 

and  fl  denote  a  floating-point  arithmetic  computation. 
Then  the  following  theorem  gives  a  bound  on  the  round¬ 
off  errors  incurred  in  performing  basic  arithmetic 


' 

. 

' 


j 


operations  in  floating-point  using  a  double-precision 
accumulator : 


16 


x 


1 


2 


b 


1 


Theorem  2.2  Given  two  floating-point  numbers 

h2 

(a1)  and  x2  =  2 


(2.17) 


fl(x2*x2)  =  (x 2*x2)(l+e)  J 


where  *  denote  any  of  the  operators  and  £ 

is  defined  by  (2.16). 


For  the  proof  of  this  theorem  see  [17,pp.7-9]  or 
[4  ,pp  .  88-89  ] . 


2.3.2  Inner  Product  Computations  Throughout  this 

paper  we  assume  that  floating-point  computations  are 

performed  with  a  precision  of  t  binary  digits  on  a 

computer  with  a  double-precision  accumulator  but  without 

accumulation  of  inner  products.  The  inner  product 
n 

l  a.b.  without  accumulation  in  floating-point  is  computed 
i=l  1  1 

as  follows :  the  products  a^b^  are  comPuked  with  2t 
digits  and  rounded  to  t  digits.  The  resulting  products 
are  then  added  in  the  order  in  which  they  are  written  using 


■ 
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a  double-precision  accumulator. 

Theorem  2 . 3  Using  a  double-precision  aooumulator 
without  floating-point  accumulation , 

n  n 


(2.18) 

i=l  1 

where 

(2.19) 

-t, 

\e2\  <  n2  , 

(2.20) 

~t  1 

lerl  f  ( n-r  +2  ) 2  ,  r=2i3i...Jni 

and 

(2.21) 

- 1  ?  - 1 

2  =  1.06  2  3 

provided 

(2.22) 

n2~tl  <0.1  . 

A  proof  of  this  theorem  is  given  in  [ 17  ,pp . 18-19 ] . 


‘ 


; 

V 


' 


18 


2.3.3  Matrix  Operations  We  will  repeatedly 
require  the  round-off  error  bounds  for  the  matrix  addi 
tion  and  multiplication.  Theorems  2.2  and  2.3  can  be 
used  to  prove  the  following  two  theorems: 

Theorem  2.4  If  A  and  B  are  matrices  of 
the  same  dimensions ,  then  using  a  double-precision 
accumu lator  3 


(2.23) 

fi(A+B)  =  A+B+F2  , 

where 


(2.24) 

\F2\  <  2~t\A+B |  . 

Theorem  2.5  If  A  and  B  are  matrices 

conformable  for  matrix  multiplication 3  then  using  a 
double-precision  accumulator  without  floating-point 
accumulation  of  inner  product 3 


fUAB)  =  AB+F2  , 


(2.25) 


' 

■■■■ 


■ 
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where 


(2.26) 

- 1 1 

\F  2  |  <  n2  d\A |  | B |  , 

and  n  is  the  common  dimension  of  A  and  B. 

The  proof  of  Theorems  2.4  and  2.5  can  also  be  found  in 
[17]. 

The  following  corollaries  are  a  direct  consequence  of 
Theorem  2.4  and  Theorem  2.5,  respectively: 


Corollary  2.4.1 

(2.27) 

IkjlL  5  • 

(2.28) 

Corollary  2.5.1 

llp2IL  ;  nS  *JIMIL  ML  • 

‘ 


l 


■■ 


■ 
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CHAPTER  III 


INVERSION  AND  SOLUTION  OP  EQUATIONS 
BY  USE  OF  PARTITIONED  MATRICES 


3 . 1  Schur's  Identity 

It  is  sometime  advisable  to  use  partitioned  matrices 
when  inverting  fairly  large  sized  matrices.  In  this 
section  we  will  consider  Schur’s  identity  [9]  for  matrix 
inversion . 

Let  R  denote  a  nonsingular  matrix  of  order  2n. 

We  write 


(3.1) 


R 


A  j  B 
C  I  D 


where  A  and  D  are  nonsingular  square  matrices  of  order 
r  and  2n-r,  respectively.  Thus  B  is  of  dimension 
r-by-(2n-r)  and  C  (2n-r)-by-r.  The  broken  lines  indicat 
matrix  partitioning.  Similarly,  we  partition  R-1  in 
exactly  the  same  manner  as  R, 

E  j  F" 

G  I  H 

Since  RR”1  =  I,  it  follows  that 


(3.2) 


R 


-1 


. 


. 


- 


■ 
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(3.3) 


f 


AE 

+ 

BG  = 

I 

r,r 

y 

AF 

+ 

BH  = 

0  0 

r , 2n-r 

y 

CE 

+ 

DG  = 

On 

2n-r , r 

y 

CF 

+ 

DH  = 

T 

V 

2n-r , 2n-r 

y 

where  the  subscripts  denote  the  orders  of  the  identity 
and  zero  matrices.  Premultiplying  the  second  equation 
of  (3.3)  by  CA-1  and  subtracting  from  the  fourth,  we 
obtain 


(3.4) 


(D-CA  1B)H 


'2n-r 


2n-r 


j 


or 


(3-5) 


3 


where 


(3.6) 


A  =  (D-CA"1B)  . 


From  the  second  equation  of  (3.3), 


.  *  hs  f 


' 

. 
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(3.7)  P  =  -A  1BH  . 

Similarly,  the  first  and  third  equations  of  (3 

(3.8)  G  =  -HCA"1  , 


and 


(3.9) 


E  =  A  1  -  A  1BG 


Thus  we  obtain 


Schur’s  Algorithm  I  To  compute  R 
R  in  (3.1), 


Step  1 : 

Compute 

(1.1) 

1 - 1 

1 

< 

(1.2) 

A  XB  , 

(1.3) 

C(A_1B)  , 

(1.4) 

A  =  D  -  C(A_1B) 

Step  2  : 

Compute  H  =  A 

Step  3  : 

Compute  P  =  -(A 

Step  4  : 

Compute 

(4.1) 

HC  , 

(4.2) 

G  =  -(HC)A"1  . 

3)  give 


from 


- 


' 


. 


v 
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Step  5 :  Compute 

(5.1)  ( A-1B) G  , 

(5.2)  E  =  A-1  -  ( A-1B)G  . 

In  practice,  all  intermediate  results  which  may  be  required 
later  are  saved.  For  example  computed  in  Step 

1.2  is  saved  to  be  used  in  Steps  1.3,  3,  and  5.1* 

3 • 2  An  Error  Analysis  of  Schur's  Algorithm  I 

In  this  section,  we  place  norm  bound  on  the  round¬ 
off  error  incurred  in  computing  R-1  using  Schur’s 
algorithm  I.  Let  ,  i=l,2,...,5  be  the  error  at  the 
i  step  and  ,  the  error  at  the  Substep  i*j  of 

Schur’s  algorithm  I.  The  subscript  e  and  c  will  be  used 
to  denote  exact  and  computed  quantity,  respectively. 

3.2.1  A  Bound  for  ||E1||oo  The  computational  equation 
for  the  calculation  of  A  in  Step  1  becomes 

(3.10)  Ac  E  [D-C{(Ag1+E11)B+E12}+E13]+Ell]  , 


or 


A  =  ( D-CA-1B )  - 
c  e 


C  E  ^  ^  B  -  C  E  ^  2  "*■  ^  3  ^  1 4 


(3.11) 


' 

■ 


, 


‘ 

. 


■ 
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But 


=  D-CA  XB  . 
e  e 


Therefore 

(3.12)  |A0-Ae|  =  |(-CE11B-CE12+E13+Ei1))  I  , 

or 

(3-13)  ll|E1|||„  =  |||CE11B+CE12-E13-E11(|||00  . 

Using  (2.5)  and  Property  (C)  in  Subsection  2.1.4,  it 
can  be  seen  that 

<(3.14)  II Ej.  <  ||C||J|E1JJ|b||„+||c||Je12||oo+||E13||oo+||ei1)||co 

We  need  satisfactory  bounds  for  the  norms  of  error  terms 
on  the  right  side  of  (3.14). 

3. 2. 1.1  A  Bound  on  ||E^1||oo  For  the  computed 

inverse  A-1  of  a  non-singular  matrix  A  of  order  n, 

0 

Wilkinson  [17]  has  shown  that 


' 


. 

, 

. 


5,1  ■  i  -  1 
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(3.15) 

> 

> 

O  I 

M 

1 

H 

II 

« 

where 


(3.16) 

IlKlL  <  eg(  2 . 005  n2+n3)||Ag1||„  , 

2 

with  terms  of  order  e  ignored.  g  denotes  the  maximum 
element -in  magnitude  at  any  stage  in  the  reduction  of  A 
into  upper-triangular  matrix.  Premultiplying  (3.15) 
by  A  ^  and  using  (3.16),  it  follows  that 


(3.17) 

II  A"1  <  eg(  2 . 005  n2+n3)||  A^H  J|  A^H  „  . 

Since  the  matrix  A  in  (3.11)  is  of  order  r,  we 
obtain  from  (3*17),  the  following  norm  bound  on  the  error 
incurred  in  the  calculation  of  inverse  of  A: 


(3.18) 

llEnIL  5  egA(2.005  r2+r3)||Ae1||J|A01||„  , 

where  gA  denotes  the  element  of  maximum  modulus  at  any 
stage  in  the  reduction  of  A  into  upper-triangular 
matrix  using  Gaussian  elimination. 


. 


' 


„ 


' 


. 
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3. 2. 1.2  A  Bound  on  ||E,J|  Let  (A  1B)  , 

_ 11  12"  «>  c  e  and 

(A  ^B)  denote  the  exact  and  computed  values  of  A-1B, 

O  C  0 

respectively . 

Then 

(3.19)  |e12I  -  Ka;1b)0  -  (A;1B)ei  . 

Using  (2.28)  and  noting  that  A-1B  is  of  dimension 

c 

r-by-(2n-r),  we  obtain  the  following  norm  bound  on  the 

error  incurred  in  the  computation  of  A_1B: 

c 

(3.20)  l|E12IL  <  erllA^lljBlI,  . 

3.2.1. 3  A  Bound  on  || E-^^ll oo  Let  (CtA^B^^ 
and  (C(A_1B)  )  denote  the  exact  and  computed  values 

C  0  0 

of  the  matrix  product  C(A”^B)  ,  respectively. 

Then 


E13l  -  I (C(A"1B)c)c  -  (C(A;1B)0)e|  . 


(3.21) 


' 


■ 


■ 
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Using  (2.28),  we  derive 


(3.22) 

l|E13IL  i  er||  C||  J  (A~1B)C»00  . 

Since 


(3.23) 

(A^B)c  =  (A;lB>e+E12  - 

it  follows  from  (3.22)  that 


(3.24) 

llEl3IL  i  er>||  C||  TO(||  A~1B||  «,+!!  E12||  „)  . 

Thus  using  (3-20)  in  (3.24),  we  obtain 


(3-25) 

l|E13IL  i  er(l+re)||c|| JA^II J|b||„  . 

seen  that 

3.2.1. 4  A  Bound  on  ||  E  ^  24 II  oo  ^  can  be  easily 
computed  value  of  CA-1B  satisfies  the  computa 

tional  equation 


(3.26) 

(CA_1B)c  =  C(Ac1B+E12)+E13  . 

(3.26) 


- 


■ 


Therefore,  it  follows  from  (2.27)  that 


(3.27)  l|ElJ4||00  <  e||D-{C( Ac1B+E12 )+E13>|| 

or 


(3.28)  l|Eli4||OT  <  £||D||eo+||C|L(||A;1||  J|B|L+||E12||J+||E13|| 
Substituting  (3.20)  and  (3.25)  in  (3-28), 


(3.29)  ||El4IL  <  £{||D||t<>+(l+re)2||C||J|A;1||  Jb^}  . 

Using  (3.18),  (3.20),  (3.25)  and  (3.29)  in  (3.11) 
we  obtain  the  following  norm  bound  on  : 

(3-30)  || Ej.  <  £||D||(o+e{gA(2.005  r2+r3 )|| A~  1|| oo+l+2r 

+r(r+2)e}||c||J|A;1||J|B||eo  , 


with  terms  of  order 


ignored . 


3.2.2  A  Bound  for  IIeJI  Let  A  ^ 

'  d1'00  0 

denote  the  exact  and  computed  values  of 


and  A-1 
c 

(D-CA"1B)'1, 


' 


- 

. 


-i 
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respectively.  Then 

(3.31)  |E.|  =  |A_1-a_1 1  . 

Also,  if  (Ag+E-L)e^  denotes  the  exact  inverse  of 
computed  (D-CA_1B),  we  may  rewrite  (3.31)  as 

(3.32)  |E2|  =  |A;1-(Ae+E1);1+(Ae+E1);1-4;1|  . 

It  immediately  follows  from  (3-32)  that 

(3.33)  IIeJL  <  Ilif-CA  +E,  jfll  +1  (A  +E.  )_1-A_1||  . 

"  2" 00  -  11  .c  e  l  e  11 00  11  e  1  e  e  11 00 

Using  (3.17)  in  (3.3*0, 

(3. 34)  ||E2||oo  <  £  gA{2.005(2n-r)2+(2n-r)3}||(Ae+E1)_1||  Ja;1!! 

+11  (A  +E-,  )-1-A-IL||  , 

where  denotes  the  element  of  maximum  modulus  in  the 

reduction  of  A  into  upper-triangular  matrix  using 
Gaussian  elimination.  The  following  theorem  is  a  result 
of  Property  (F)  in  Subsection  2.1.4  and  (3.34): 


' 


. 
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Theorem  3 . 1  Let  E denote  the  matrix  of 
round-off  errors  in  the  calculation  of  H  using  Schur’s 
algorithm  I  (see  Section  3.1),  then  if  II  A~  ^||  Jl  S-jll  „  <  1, 


(3.35)  Us 


< 

00  — 


■  g  ^{2 . 0  05 ( 2n-r )2  + ( 2n-r ) 3 


}||a:jii 


00 


HIa^II  JSjII 


where  ||£’1||  is  defined  by  (3.  SO). 
1  00 


3.2.3  A  Bound  for  ||  E  ^  ^  The  computational  equation 
for  the  calculation  of  F  becomes 


(3.36)  Pc  E  -  [{(A"1+E11)B+E12}(He+E2)+E3:L]  , 

where  denotes  the  error  in  matrix  multiplication 

(A”^B)  H  .  Since  A”1  =  He+E2 3  we  obtain  f^0111  (3.36), 

(3-37)  | E3 |  =  |Ag1BE2+(E11B+E12)A‘1+E31|  , 

or 


(3.38)  ||  E  J  <  ||  A"1)!  J|B||  J|E21l„ 


+  (I|e11||J|b||co+||e12|L)||a;1|| 


12"  00 


■c1' 


+ 


Eon|| 
31" 00 


■ 

■ 


A- 


- 


. 
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3. 2. 3.1  A  Bound  on  ||E0J|  Let  ((A  1B)  A-1 ) ^ 

ii  3 1 1  00  c  cc  e 

and  ((A  ^B)  A"1)  denote  the  exact  and  computed  values 
0  c  c  c 

of  (A J'B)  A  \  respectively.  Then 


(3.39) 


|E 


31 


=  I ( (A  1B)  A-1)  -  ( (A-1B)  A  1) 

1  c  cc  c  c  c  c  '  e 


Using  (2.28), 


(3.40)  ||E31||oo  <  E(2n-r)(||Ac1|| Jb||oo+||E12||o5)||&;1||[o  . 
Substituting  (3.20)  in  (3.40), 

(3.41)  ||E31||oo  <  e(2n-r)(l+re)||A”1||co||B||  J|  A^1^  . 

Using  (3.18),  (3.20)  and  (3.41)  in  (3.38), 

(3.42)  ||E3||00  <  llA^H  jB||J|E2||oo+e{gA(2.005  rV)^;1 


+  2n+r  (  2n-r )  e}||  Aq  1||  J|  B||  J  A"1! 


Hence,  we  state 


' 
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Theorem  3.2  Let  E denote  the  matrix  of  round¬ 
off  errors  in  the  calculation  of  F  using  Schur's  Algorithm 
I  (see  Section  3.1)3  then 

(3. 43)  ||EsIL  <  UfWjBWjB^^eig^Z.OOS  r hr 3>  ||A'  3||  „ 

+  2n+r(2n-r)z)\\A~o1\\  J|s||  Ja^IL  , 

where  II^Hoo  ^ 8  defined  by  (3.35). 

3.2.4  A  Bound  for  ||  ||  ^  The  computational  equation 

for  the  calculation  of  G  is 

(3.44)  G0  i  -  [{(He+E2)C+E1)1}(A”1+E11)+E42] 

Since  H  =  A-'*'  and  aT  5  A^^+E..,  therefore 
e  e  c  ell 

(3.45)  I E4  I  =  |A"1CE11+(E2C+Ei)1)A"1+E42l  , 
or 


I|E4IL  i  IIa;1IIJc||Je11IL+(||e2||Jc||oo+i|e4i||oo)||a;1||oo+i|e42|| 


(3.46) 


' 


* 


■ 


.  • 

. 


■ 


, 
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3. 2. 4.1  Bounds  on  and  ||  E  ^  2 II  oo  A  norm 

bound  on  the  error  incurred  in  the  matrix  multiplication 
HcC  can  be  derived  using  (2.28).  Thus 


(3.47) 

llElJL  i  £  ( 2n-r )  ||  A c  'L|(  J|  C||  —  . 

Similarly,  a  norm  bound  on  the  error  incurred  in  the  matrix 
multiplication  (H  C)A-1  follows  from  Subsection  3-2. 1.3* 


Thus 

(3.48) 

IIe42IL  i  ^(11^41  JcIL+llEltllL)||A;1||o)  . 

Using 

(3.47), 

(3.49) 

lEH2i«  ^  er (1+  ( 2n-r ) e > ||  A~ 1 1| J  C||  J  A”4l »  . 

The  following  theorem  is  a  result  of  (3.18),  (3.46), 


(3.47) 

and  (3.49)  : 

Theorem  3.3  Let  E A  denote  the  matrix  of  round- 

off  errors  in  tine  calculation  of  G  using  Schur’s 
Algorithm  I  (see  Section  3.1 ),  then 


■ 


■ 


J1  i!  !  Ut  ■'  '  ■  «3ll  ■  '  -;) 


' 

' 


' 


- 
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(3-50)  IlffJL  <  i\\E2\'+t{gA(2.00l>  r‘2+v3)\A-e1\\Jh-e1\\J 

+  e{2n+r  ( 2n-r)  e}\\  A ~  J |j  ^ ] ||  j4 ~  J  j|  J  C||  —  , 

where  ||  «>  is  given  by  (3.35). 

3.2.5  A  Bound  on  || E ^ ^  The  computational  equation  for 
the  calculation  of  E  is 


(3.51) 

or 

(3.52) 

or 

(3.53) 


5  [a;1+eh 


-  {(Ae1+E11)B+E12}(Ge+E1))+E51+E52 


Ell-AelBE4+(EllB+E12)GC+E5l+E52 


» 


l|E5L  <  iBnL+lA^ljBljEnH.+  dE^ljBl.+  lE^I.JlO, 


+  I|e51IL+||e52|| 


We  now  proceed  to  place  norm  bounds  on  the  matrices 
and  Eco. 


3. 2. 5.1  Bounds  on  llE^-JI^  and  ||E^2II00  Norm  bounds 


for  the  matrices  of  round-off  errors  E._,  and  E,-0  incurred 

51  52 

in  the  calculation  of  (A_1B)G  and  A-1  -  (A-1B)G,  respectively 


- 


>:  '  ■  r  '  >  •  f  : 
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can  be  obtained  by  using  an  analysis  analogous  to  that  in 
Subsection  3-2.1.  We  derive 

(3.54)  ii e51II oo  ;  e(2n-r)(l+re)||A"1||  Jb||JGc||oo  , 
and 

(3-55)  ||E52||oo  <  £[l+(l+re){l+(2n-r)£}||B||oo||GcIL]||A;1||t<>  . 

The  following  theorem  is  a  result  of  (3-18),  (3-20), 
(3-53),  (3-54)  and  (3.55)*: 

Theorem  3-^  Let  E ^  denote  the  matrix  of 
round-off  errors  in  the  oaloulation  of  E  using  Sohur's 
algorithm  I  (see  Section  3.1),  then  if  terms  of  order 
e2  are  neglected, 

(3-56)  ||ZSIL  <  \\A-e1\\jB\\jE4\\m+z{l+gA(2.005  r2 +r*  )  ||  4' 2|| 
x  WbfWj-ztg/Z-  005  r2+rS  )\\A'e1\\m+l  +  2n(l  +  e) 

+  rtZn-rl^WA-fW  J|  S||  J  ffJL  , 

where 


is  defined  by  (3.50). 


- 

M  I  ^  mm  I 


. 


' 


’ 

. 


£ '  I 
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3 • 3  A  Bound  for  the  Computed  Inverse  of  R 

-1  -1 

Let  R  and  R  denote  the  exact  and  computed 

G  0 

inverse  of  R,  respectively.  Then  the  computational  equation 
for  the'  calculation  of  R-"*"  becomes 


(3.57) 


E  +E,-  F  +E0 
e  5  e  3 


G  +E i.  i  H  +E0 
e  4  i  e  2 


3 


where  E^,  E^,  E^  and  E,_  are  given  by  (3*35),  (3.43) 

(3.50)  and  (3.56),  respectively. 

E  I  F 
e  1  e 

__  -1 

Since  R0  =  - j -  , 


it  follows  from 

(3.58) 


(3.57)  that 


Hence , 


(3.59)  || r01  -  Rg1!! i  IIe2IL+IIe3IL+IIeiiL+IIe5IL  . 

where  HE^,  ||E3||ra,  |1  E^H „  and  ||E5||m  are  bounded  by  (3-35), 
(3.43),  (3.50)  and  (3-56),  respectively. 


' 


* 


' 

2a  - 
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3.^  A  Smaller  Bound  for  II R  -  R  "Hi 
_ 11  c _ e  11  °° 

A  smaller  bound  on  ||  R~ 1  -  R”‘*'||oc)  can  be  placed  by 
using  the  definition  of  “-norm.  Since  the  maximum  row- 
sum  of  error  matrix  on  the  right  side  of  (3*58)  can 
occur  either  in  matrix  [E^ |  E^]  or  [E^ j  E^],  it 
immediately  follows  that  at  least  one  of  the  following 
is  true: 


(3*60) 


3 


This  bound  is  at  least  as  small  as  that  given  by  (3.59). 


3 . 5  Solution  of  System  of  Equations 
Let 


(3.61) 


Rx  =  b  , 


be  a  system  of  linear  algebraic  equations  of  order  2n 
with  nonsingular  coefficient  matrix  R.  Then,  we  can 
represent  (3.61)  in  the  form 


"A 

_X1 

~bi 

C 

D 

_X2_ 

1  O'  1 

rv>  l 

1 - 

(3.62) 


3 


' 

■ 


«  *  ;  -  ■ 


. 


' 

. 

1 

■ 


■ 
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where  A  and  D  are  nonsingular  square  matrices  of 
order  r  and  2n-r,  respectively.  The  matrix  B  is  of 
dimension  r-by-(2n-r)  and  C  of  dimension  (2n-r)-by-r. 
Thus  the  vectors  and  b^  are  of  order  r  and  vectors 

x2  and  b^  of  order  2n-r,  each.  The  broken  lines 
indicate  matrix  partitioning.  Equating  terms  on  both  sides 
of  (3.62), 

r 

*  Ax.^  +  Bx2  =  b^  , 

(3.63)  \ 

I  Cx1  +  Dx2  =  b2  . 

Eliminating  x1  from  (3.63), 

(3.64)  (D-CA_1B)x2  =  b2-CA"1b1  . 

Thus,  if  the  matrix  (D-CA_1B)  is  nonsingular,  the  system 
(3.64)  can  be  solved  for  x2  using  Gaussian  elimination. 
The  solution  vector  x2  can  not,  however,  be  substituted 
in  the  second  equation  to  determine  x-^,  since  the  matrix 
C  is  rectangular  for  all  r  t  n.  This  suggests  that,  in 
such  case,  we  resort  to  a  method  of  evolving  a  system  of 
equations  similar  to  (3.64).  Thus  eliminating  x2  from 
(3.63)  , 


. 


' 

' 


■ 


' 

. 


.. 
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(3.65)  (A-BD  1C)x1  =  b1-BD  Xb2  . 

Hence,  we  obtain 

Schur's  Algorithm  II  To  compute  x=(x1,x2) 
from  (3.63)  when  r  ^  n 


Step  1: 

Compute 

(1.1) 

A'1  . 

(1.2) 

CA"'1  , 

(1.3) 

( CA-1 ) B  . 

Step  2 : 

Compute 

(2.1) 

(CA"1)b1  , 

(2.2) 

bp-(CA-1)b,  . 

Step  3: 

Compute  x2  from  (3.64). 

Step  4 : 

Compute 

(4.1) 

D"1  , 

(4.2) 

BD"1  , 

(4.3) 

( BD-1 ) C  . 

(4.4) 

A  -  ( BD-1 ) C  . 

Step  5: 

Compute 

(5.1) 

(BD_1)bp  , 

(5.2)  b1-(BD  1 )b 2  . 

Step  6:  Compute  x-^  from  (3.65). 


. 


- 


' 

V 


■ 

' 


When  the  matrices  A,  B,  C  and  D  are  square,  the  last 
three  steps  in  the  above  algorithm  are  replaced  by 
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Step  4 :  Compute 


(4.1)  Dx2  , 

(4.2)  b^-Dx^  . 


Step  5:  Compute  x-^  from  Cx-^  =  b2~Dx2  . 


3 • 6  An  Error  Analysis  of  Schur's  Algorithm  II  with  r  ^  n 

We  now  proceed  to  place  norm  bounds  on  the  round-off 
error  incurred  in  computing  the  solution  of  the  system  of 
equations  (3.61)  using  Schur's  Algorithm  II.  We  will 
use  the  same  notations  as  in  Section  3*2. 

3.6.1  Bounds  for  ||E1||oo  and  ||E2||oo  A  norm  bound 
for  the  error  in  computing  D-CA  '*'B  has  already  been 
obtained  in  Subsection  3.2.1.  Although  the  algorithm  for 
calculating  D-CA-1B  in  Section  3.1  is  slightly  different 
from  that  in  Section  3.6,  it  can  be  easily  shown  that  round¬ 
off  error  bounds  are  exactly  the  same.  Therefore, 


(3.66) 


llElL  1  e||D||<<1+£{gA(2.005  r2+r3)||Ae1||„+l+2r 


ope 


' 


■.  ii 

-  —  —  - 


■ 


i  3-  »  •  e  l[s  >  ix 
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with  terms 
Similarly , 
(3.66)  by 


3 

of  order  e  ignored. 

a  norm  bound  for  can  be  derived  from 

replacing  D  by  b^,  and  B  by  b^  , 


(3.67) 


e2IL  <  £l|b2l|oo+e{gA(2.005  r2  +  r3)||  Ae1||oo+l+2r 


+ 


r(r+2)e}||  C 


A 


-1 


3.6.2  A  Bound  for  ||  E  ^  ^  Wilkinson  [17]  has  shown 
that  the  computed  solution  of  Ax=b ,  is  an  exact  solution 
of  the  perturbed  system  (A+K)x=b, 


where 

(3.68)  || K|| m  <  eg( 2 . 005  n2+n3)  . 


g  denotes  the  element  of  maximum  modulus  at  all  stages 

in  the  reduction  of  A  using  Gaussian  elimination.  If 
(2)  (2) 

x  and  x  denote  the  exact  and  computed  solution 

e  c 

(2) 

of  (3.64),  it  follows  that  x  is  an  exact  solution 

c 

of  the  perturbed  system 


<VK1)xo2)  =  be2)+E2 


(3.69) 


y 


' 


■ 


-  n  !_  f  I  b  n  ■  : 


. 


and 


where  b^2^  denotes  the  exact  value  of  b 
(3-70)  IlK-Jco  <  egA{2.005(2n-r)2+(2n-r)^}  . 
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gA  denotes  the  element  of  maximum  magnitude  at  all  stages 

in  the  reduction  of  A  into  upper-triangular  matrix  using 

Gaussian  elimination.  Since  A  =A  +E. 

c  e 

we  obtain  from  (3*69), 


j-.  and  A  x^2^=b^2^ 
1  e  e  e 


(3.71) 


(2)-xf2)|  =  |Ao1E,-a=1(E1+K1)x 


(2) 


or 


(3.72)  ||e j  <  «e2|| Ja;1il+(||e1||oo+||k1||ji|a;1|Jx(2)|| 


since  | E  |  =  |x(2)-x(2) | 

1  3  1  c  e  1  . 

3.6.3  Bounds  on  ||  E  ^  ^  ,  ||  E^||  ^  and  ||  Eg  ||  ^  A  reference 
to  equations  (3.66),  (3.67),  (3.70)  and  (3.72)  shows 
that  norm  bounds  for  E^ ,  Eg  and  Eg  can  be  obtained 
merely  by  replacing  D  by  A,  C  by  B,  A  by  D,  B  by  C, 

A”1  by  6”1,  b2  by  b1?  by  b2,  (2n-r)  by  r,  r  by  (2n-r), 
x £2^  by  x^,  E±  by  E^ ,  E2  by  Eg ,  Eg  by  Eg  and  K±  by  K2  . 


. 


■ 


' 


* 


. 


■ 


■ 
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Therefore , 

(3-73)  ll e6 | „  <  |e5ij«;1|.+(|e1,|„+|k2|ji«;1|Jx(1)||1„  , 

where 

(3.74)  llE^H^  <  e||A||oo+e[gD{2.005(2n-r)2+(2n-r)3}||D^1||oo 

+  1+2  (  2n-r )  +  (  2n-r )  (  2n-r+2 )  e  ]||  C||  J  D~ X||  J|  B||  w  , 

(3.75)  Halloo  <  £||b1||oo+e'[gD{2.005(2n-r)2+(2n-r)3}||D”1||oo 

+  1+2  (  2n-r )  +  ( 2n-r )  (  2n-r+2 )  e  ]  ||  B||  J|  D"^  J  b2||  „  , 

and 

(3-76)  || K2|| <  e  g5(2.005  r2+r3)  . 

6  and  denote  the  exact  value  of  A-BD~2C  and  the 

e  c 

computed  solution  of  (3.65),  respectively. 


The  following  theorem  is  a  result  of  (3.72)  and  (3.73) 


' 


■r  .  1  '• 


■ 


■ 


. 


■ 


1 
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Theorem  3-5  Let  x  and  x  denote  the  exact 

- - —  e  c 

and  computed  solution  of  (3.61)  using  Schur's  Algorithm  II 
(see  Section  3.5),  then  at  least  one  of  the  following  is 
true  : 


r 


-It  n  (2) 


(3.77)  |ka-*eL  5  < 


II  ss« «»  3 « (I  *4 1  „+l  X2I J I  J I  J  ' «  »  - 


-III  II  (1) 


V 


where  norm  hounds  on  matrices  E  ,  E  ,  E  A,  E K1  and 

1  <u  4  O  1 

Kg  are  given  hy  (3.66),  (3.67),  (3.74),  (3.75),  (3.70)  and 

(3.76),  respectively. 


3 . 7  An  Error  Analysis  of  Schur’s  Algorithm  II  with  r=n 
In  the  following  we  propose  to  place  norm  bounds  on 
the  matrices  of  round-off  errors  and  incurred 

in  computing  the  last  two  steps  of  Schur's  Algorithm  II 
with  r=n  (see  Section  3.5).  The  error  analysis  of  the  first 
three  steps  of  this  algorithm  has  already  been  given  in 
Section  3.6. 


3.7.1  A  Bound  on  ||  E  ^  ||  ^  The  computational  equation 
for  the  calculation  of  b2~Dx2  is 

=  {b2-Dx£  ^+Ei|1>+E2|2  , 


(3.78) 


b  2~Dx2 


' 


.  .  -  *  '•  c 


■ 


. 


r 

• 

' 


' 


'Wtf 
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or 

(3.79)  b0-Dx0  =  {b0-D(x^2^+x^2^-x^2^  )+E,n  }+Ejlo  , 

d  d  d  e  c  e  41  h  d  ’ 

from  which  it  follows  that 

(3.80)  |E4IL  <  IlDllJx^^x^^lL+llE^IL+llE^L  . 

We  now  need  satisfactory  norm  bounds  for  the  matrices 

and  E ^  2  .  A  bound  for  ||  x^ 2  ^  -x^ 2  ^  ||  ^  is  given  by 

(3-72),  since  E0  =  x^^ -x^^  . 

3  3  c  e 

3. 7. 1.1  Bounds  for  [|E1|1||oo  and  ||E1|2||oo 
Using  (2.27)  and  (2.28), 

(3.81)  |e4iIL  <  £  nil D||  , 

and 

(3.82)  ||E42S„  <  e||b2-Dx(2)+E4l|L  • 

Using  (3-81)  in  (3*82), 

l|E42L  i  £{||b2!|„+(l+ne)||D||  Jx<2)fl_)  . 


(3.83) 


■ 


' 


* 


' 
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Substituting  (3.81)  and  (3-83)  into  (3*80),  we 
obtain  the  following  norm  bound  for  the  matrix  of  round¬ 
off  errors  incurred  in  computing  b^-Dx^: 

(3.84)  || <  e[||b2||co+(n+l+n£)||D||co!|x^2)||0O]+l|D|LI|x^2)-x^2)|l 

3.7-2  A  Bound  on  || E^||  ^  Using  a  method  analogous  to 
that  in  Subsection  3.6.2,  we  derive  that 


where 


(3.86)  |  <  e  gc(2.005  n2+n3)  . 

A  bound  for  JeJI^  is  given  by  (3.84). 

The  following  theorem  is  a  result  of  (3-72)  and 

(3.85) : 


Theorem  3.6  Let  x  and  x  denote  the  exact 

- — —  g  c 

and  computed  solution  of  (3.61),  then  if  the  matrices 
A,  B,  C ,  D  are  square,  at  least  one  of  the  following  is 


■ 


II  s  ci|,  i|a :|+ :  (  •  „:a||(3nv>n)+„;;:  /  n:  < Pti •  ’ 


- - — i  — -  • 


■ 


. 


. 
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true 


(3.87) 


£■  „||  ||&”il|  +(!|e,||  +|| a:, II  )||a"j||  ||x  „ 

2"  00*1  Q  ‘I  00  •  ]l'00  I!  2  00  00 1'  (2  *'00 


-III  II  (2) 


x  -x  < 

a  e 11  oo  _ 


^  )  ||  I  ||  Q~  1 1 

^11  00  '  II  2 ^  00^  **  C  00  ll  oo  * 


oi^^ii  oo+ii-^oii  j*; 


where  norm  hounds  on  E  ^3  E^3  E^3  Kj  and  K ^  are  given  by 
(3.66)3  (3.67)3  (3.74)3  (3.70)  and  (3.76)3  respectively 3 

with  r  replaced  by  n. 


3 . 8  Comparison  of  Norm  Bounds 

A  reasonable  theoretical  comparison  of  the  norm  bounds 
on  the  round-off  errors  incurred,  in  computing  the  inverse 
of  a  matrix  and  solution  of  a  system  of  linear  algebraic 
equations,  using  Gaussian  elimination  and  Schur’s  Algorithms 
(see  Section  3.1  and  Section  3.5),  can  not  be  done  due  to 
the  nature  of  terms  involved  in  the  latter.  However,  a 
comparison  of  the  norm  bounds  when  r=n,  can  be  made. 

Without  loss  of  generality,  we  can  assume  that  elements 
of  the  matrix  R  have  been  scaled  such  that 


(3.88) 


Therefore , 


' 


. .  „  ..  -  ;  ■  ^ 


- 


■ 


■ 
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(3-89)  ilAlLJBjl^JCll^  and  i|  D||  _  <  n  . 

From  the  definition  of  °°-norm  of  a  matrix  and  (3.57) 
it  follows  that 


(3.90) 

K1IL  £  IIR^IL  » 

(3.9D 

iia;xl  s  iir;1il  . 

(3.92) 

I|Q0L  £  IIR^IL  . 

(3-93) 

I|ECIL  £  !Iro1IL  » 

and 

(3.94) 

ii  Eeiioo  £  iir;xil  • 

Also,  since 

E  =  A"1-(A_1B)G, 

(3.95) 

IIa^IL  £  Hr^IL  . 

and 

(3.96) 

IIA^IL  £  ! i  R q  1 II  co  • 

' 


‘ 


- 
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Using  inequalities  (3.89)  to  (3.96)  in  (3*30),  (3-35), 
(3.42),  (3.50),  (3-56)  and  (3.60),  it  follows  that  the 
dominant  term  in  the  expression  for  round-off  error  incurred 
in  computing  the  inverse  using  Schur's  algorithm  I  (see 
Section  3.1)  is  0(n  )  a  >  5,  whereas  the  corresponding 
error  term  using  Gaussian  elimination  is  0(n  ).  The  latter 
follows  from  (3.17)  by  replacing  A  by  R  and  n  by  2n. 
Similarly,  it  can  be  shown  that  the  dominant  term  in  the 
expression  for  round-off  error  incurred  in  computing  the 
solution  of  a  system  of  linear  equations  using  Schur’s 
algorithm  II  (see  Section  3*5)  and  Gaussian  elimination  is 
0(n  ),  a  >  5,  and  0(n  ),  respectively.  However,  the  actual 
errors  depend  on  the  nature  of  the  elements  of  the  matrices 
considered . 
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CHAPTER  IV 


BLOCK-SYMMETRIC  MATRICES 

4 . 1  Inversion  by  Use  of  Schur's  Identity 

Let  R  be  a  nonsingular  block-symmetric  matrix  of 
order  2n.  We  write 


(4.1) 


R  = 


A  B 


B  I  A 


where  A  and  B  are  square  -matrices  of  order  n. 
Similarly  we  partition  R-1  in  exactly  the  same  manner 
as  R, 


(4.2) 


R 


-1 


E  ]  F 
F  !  E 


It  can  be  easily  shown  by  equating  the  terms  on  both  sides 
of  RR-1  =  I,  that 


(4.3) 


E  =  (A-BA  1B)  1  , 


(4.4)  F  =  (B-AB  1A)  1  , 

provided  the  matrices  A,  B,  (A-BA  ^B)  and  (B-AB  "''A)  are 
nonsingular . 


Thus  we  obtain 


' 


.  ■» 


' 
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Schur's  Algorithm  III  To  compute  R  1  from 
R  in  (4.1), 


Step  1 : 

Compute 

(1.1) 

A'1  . 

(1.2) 

A_1B  , 

(1.3) 

B ( A_1B )  , 

(1.4) 

A  =  A-B ( A 

Step  2  : 

Compute 

Step  3  : 

Compute 

(3.1) 

__  #\ 

i — 1 

i 

m 

(3.2) 

B_1A  , 

(3.3) 

A(B"1A)  , 

( 3 . 4) 

6  =  B-A ( B 

Step  4: 

Compute 

4 . 2  An  Error  Analysis  of  Schur's  Algorithm  III 

In  this  section  we  place  norm  bounds  on  the  round¬ 
off  errors  incurred  at  each  step  in  the  computation  of 
R  ^  using  Schur’s  Algorithm  III.  We  will  use  the  same 
notations  as  in  Section  3.2. 


4.2.1  Bounds  on  ||E1||oo,  ||E2||oo,  II  E^  and  He^ 
The  error  analysis  of  the  inversion  of  a  block-symmetric 
matrix  using  Schur's  Algorithm  III  is  similar  to  that 


‘ 


■ 


' 
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developed  in  Section  3.2.  Bounds  for  ||  E-.  ||  ^  and  ||  E  ^  ^ 
follow  from  (3-30)  and  (3-35),  respectively,  by 
replacing  r  by  n,  C  by  B  and  D  by  A.  We  have 
therefore 


(4.5) 

IIEJL  i  £||A||oo+£{gA(2.005  n2+n3)||Ae1||oo+l+2n 

+  n(n+2)e}||B||  Ja;1!  JBIL  , 

•5 

with  terms  of  order  z  ignored,  and 


(4.6) 

{Eg  (2.005  n^+n3 )  ||  A  1|  +||  A  1||  J|  E^||  }  ||  A  1|| 

IIe2IL  i  — - 2 - 5 — - — 5 — 

i-ll  a;1!  JeJL 

Bounds 

for  ||  E^||  oo  and  ||  E^||  ^  follow  similarly  from 

(3.30) 

and  (3*35) a  respectively,  by  replacing  r  by  n. 

A  by  B,  B  by  A,  C  by  A,  D  by  B,  A  by  6  and 


E1  by 

E^,  hence. 

(4.7) 

||E3IL  <  e|| B||<io+e{gB(2 .005  n2+n3)||B;1||oo+l+2n+n(n+2)e} 

x  II a|| Jb;1!!  jail  , 

with  terms  of  order  z  ignored,  and 


■  ■ 


' 


■ 


. 


- 
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{egg(2.005  n2+n^)||5c1 


00 


(4.8) 


< 


1- 


3" 00 


The  following  theorem  is  a  result  of  (4.1),  (4.2),  (4.6) 
and  (4.8)  : 


Theorem  4.1  Let  R~ 1  and  R~  1  denote  the 


e  o 


exact  and  computed  inverse  of  a  nonsingular  block- symmetric 


A  i  B 

matrix  R  =  of  order  2n3  where  A  and  B  are 

square  matrices  of  order  n.  If  the  matrices  A ,  B, 

-1  - 1 

A-BA  B  and  B-AB  A  are  nonsingular ,  then  the  norm 


bound  for  the  round-off  errors  incurred  in  computing  R 
using  Schur’s  Algorithm  III  (see  Section  4.1)  is  given  by 


(4.9) 


where  Halloo  and  ||£'^||oo  are  given  by  (4.6)  and  (4.8) , 
respective ly . 

4 . 3  Inversion  by  Use  of  Charmonman's  Algorithm 

A  method  for  computing  the  inverse  of  a  block-symmetric 
matrix  has  been  proposed  by  Charmonman  [1].  This  method 
is  more  efficient  in  the  sense  that  it  results  in  saving 


. 


— 


. 

■ 


A 


, 

' 
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of  computer  storage  and  the  number  of  arithmetic  operations. 
The  algorithm  given  by  Charmonman  [1]  can  be  devised  by 
use  of  the  following  theorem: 


A  i  B 

E  !  F 

Theorem  4.2  The  inverse  of  P  = 

_  _i_  _ 

B  1  A 

is 

I 

i 

i 

i 

,  C*j  I 

where 


(4.10) 


r 

E  =  0 . 5 (P~1+Q~1 )  } 
F  =  0 . 5 ( P~ 1 -Q~ 1 )  , 

P  =  A+B  , 


Q  =  A-B  . 

V 


Proof :  The  proof  of  the  first  part  of  the 

theorem  follows  immediately  from  Schur’s  Algorithm  I 
Section  3.1)  by  replacing  C  by  B  and  D  by  A. 
second  part  of  the  theorem  can  be  proved  as  follows : 
We  may  write  R~^R  =  I  as 


(4.11) 


A 

B 


B 

A 


E 

P 


F 

E 


=  I 


2n,2n  5 


(see 

The 


' 
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where  I„  _ 

2n ,  2n 

Equating  terms 


is  an  identity  matrix  of  order  2n. 
on  both  sides  of  (4.11),  we  obtain 


(4.12) 


AE  +  BF  =  I 

n ,  n 


and 


(4.13)  BE  +  AF  =  0 

n  ,n 

Adding  (4.12)  and  (4.13),  and  premultiplying  by  (A+B)-'*' 
on  both  sides,  we  derive 

(4.14)  E+P  =  (A+B)"1  . 


Similarly,  subtracting  (4.13)  from  (4.12),  and 
premultiplying  by  (A-B)-1  on  both  sides,  we  obtain 


(4.15) 


E-P  =  (A-B)  1 


The  result  in  (4.10)  follows  by  solving  (4.14)  and 

(4.15)  for  E  and  F. 


Q  .  E  .  D  . 


• 
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Thus  we  obtain 


Charmonman 1 s  Algorithm  I  To  compute  R  ^  from 
R  in  (4.1), 


Step 

1 : 

Compute 

P  = 

A+B  . 

Step 

2  : 

Compute 

Q  = 

P-2B  . 

Step 

3: 

Compute 

P_1 

• 

Step 

4: 

Compute 

Q-1 

• 

Step 

5: 

Compute 

E  = 

0 . 5 ( P-1+Q-1 )  . 

Step 

6: 

Compute 

F  = 

E-Q"1  . 

Multiplication  by  2 

and 

0 . 5  are 

done 

by  shifting  without 

any 

round-off  error. 

4.4 

An  Error  Analysis 

of  Charmonman 

’  s  Algorithm  I 

We  now  proceed 

to 

place  upper  norm  bound  on  the  round 

1  . 

off 

error  incurred 

in 

computing 

R 

using  the  above 

algorithm.  We  shall  again  use  the  same  notations  as  in 
Section  3.2. 


4.4.1  Bounds  on  IjE-JI^  and  ||  E  2  ]|  ^  Let  and  P' 

denote  the  exact  and  computed  value  of  P,  then 


(4.16) 


‘ 


S 


. 
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It  follows  from  (2.27)  that 

("•IT)  IIeJL  <  edlAlL+llBlI.)  . 

Similarly,  if  Q"  and  Q"  denote  the  exact  and  computed 

G  0 

values  of  Q, 

(4.18)  | E2 |  =  |Q'-Q'|  . 

From  the  computational  equation  for  the  calculation  of  Q, 

(4.19)  Q'  =  P^+E;l-2B+E21  , 

and  (4.18),  we  have 

(4.20)  I|e2IL  =  !|e1+e21||„  . 

or 

(4.21)  l|E2IL  <  llE1iL+l|E21l|„  . 


■ 


■ 


' 


. 


■ 
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4. 4. 1.1  A  Bound  on  ]|  ^ 


Since  E 


21 


is 


the  matrix  of  round-off  errors  incurred  due  to  the 
addition  of  computed  P  and  -2B,  it  follows  from  (2.27) 


that  ||E21||eo  <  e||p'+E1-2B||oo  .  Using  (4.17)  , 


(4.22) 


Substituting  (4.17)  and  (4.22)  into  (4.21)  give 


(4.23) 


2 

with  terms  of  order  e  ignored. 


4.4.2  Bounds  on  jj  E  3  !|  ^  an<^  ;j  E  ^  !j  ^  A  reference  to 
the  error  analysis  in  Subsection  3.2.2  shows  that  bounds 
for  !|E3||oo  and  || E ^ ^  can  be  obtained  from  (3.35) 
by  replacing  r  by  n,  A  by  P  and  Q,  respectively. 
Thus 


' 


- 


■ 

' 


-  1 
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and 

e  {gQ(  2 . 005  n2+n3 )  ||  Q^ll  „+!  Q^H  J  Ej  J  ||  Q"1!! 
(1»-25)  l|E„|L  <  - 

1  -  iiq;1iije2il 


provided 

C1*- 26)  IP^IJBJ.  <  i  . 

and 

('*•27)  iiq;1iije2il  <  i  . 

Norm  bounds  on  E-^  and  are  given  by  (4.17)  and 

(4.23),  respectively. 

4.4.3  A  Bound  on  j|  E  ^  j|  ^  The  computational  equation 
for  the  calculation  of  E  is 

(4.28)  Ec  =  0.5((Pg1+E3)  +  (Q^1+Ei|)}+E51  . 

Therefore 


(4.29) 


»e5IL  i  o-5(||e3L+|e1|II„)+||e51|| 


- 


- 

.  I  1 ■  H  1 1 "  •  > ,  - 


■ 


■ 
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4.4. 3-1  A  Bound  on  ||  E  ^  J|  ^  Since  is  the 

matrix  of  round-off  errors  incurred  due  to  the  addition  of 
P  1  and  Q  and  multiplication  by  .5  is  done  by 

V-/  O 

shifting,  it  follows  from  (2.27)  that 

('*•30)  l|E51IL  <0.5  edlP^IL+HQ^IL)  • 

The  following  theorem  is  a  result  of  (4.17), 
(4.23),  (4.24),  (4.25),  (4.26)  and  (4.27): 


Theorem  4 . 3  Let  E denote  the  matrix  of  round¬ 
off  errors  in  the  calculation  of  E  using  Charmonman's 
Algorithm  I  (see  Section  4.3),  then 


(4.3D 


g  (2.00b  n2+n3 

Ilfi'JL  <  h  e 

1+  -£ - % - - 

*  5 1  00  - 

l  2-edML+NUIlp;  !L  J 

II  Ol 


g  (2.00b  n+n^WQ-JW^  "] 

— - - — j—  HI 

1- 2z(\\A\\  „+2||s||J||  Qf  ||  J 

dMii^iigiLdip;3uj  p;2il 

2- edUL+8BIUII?e2L 


' 


■ 


■ 
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2(jMiL*2MU!lQl  JlQl 

i-2<i(\\A\\a+2\\  b||j||  of  IL 


provided 

('*•32)  ^IUIL+I|S||J||P^IL  <  2  , 

and 

('••33)  ZeClMIL+^llfilUlle^lL  <  2  • 

4.4.4  A  Bound  on  jlE^jj^  The  computational  equation 
for  the  calculation  of  F  is 

(4. 3*0  Pc  =  E+E5-Q"1+E1)  . 

Therefore 

(4-35)  l|E6IL  <  ||e4|L+||e5||„  . 

The  following  theorem  is  a  result  of  (4.25)  and  (4.31): 

Theorem  4 . 4  If  E denotes  the  matrix  of  round¬ 
off  errors  incurred  in  the  calculation  of  F  using 
Charmonman’s  Algorithm  I  (see  Section  4.3)  subject  to  (4.32) 
and  ( 4 .33) 3  then 


' 


V 


■ 


■ 

.  \ 
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(4.36) 


llE<jlL  i  \\E5\\m+z{gQ(2-00S  nhn^WQ-^W^ 

_  l 

+  2f|MIL+2||S||raJ}  _ ^  8g  _ 

2-2e<'IUIL+2||s||oo;||«^|| 


4.4.5  A  Bound  for  Computed  Inverse  of  R  If  R 
,-l 


-1 


and  R  denote  the  exact  and  computed  inverse  of  a  non- 

°  fA  !  S’ 

singular  block-symmetric  matrix  R  =  — | — 


L 


B  !  A 


of  order  2n, 


where  A 
A+B  and 

(4.37) 


and  B  are  of  order  n  each,  and  if  matrices 
A-B  are  nonsingular,  we  may  write 

E  *  F~ 

i 

I  > 

F  !  E 


and 


(4.38) 


E+Ej- 

5 

P+E6~ 

F+E6  1 

E  +  E,- 
5 

where  Er  and  Ec  denote 

5  6 

incurred  in  the  calculation 
From  (4.37)  and  (4.38) 

(*t.39)  II  r;1-^1!! 


the  matrices  of 
of  E  and  F, 
it  follows  that 


round-off  errors 
respectively . 


> 


where  E,_  and  E^  are  bounded  by  (4.31)  and  (4.36) 
subject  to  the  conditions  (4.32)  and  (4.33). 


. 


, 

. 


, 


I '  m , 


. 
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4 . 5  Comparison  of  Upper  Norm  Bounds  in  Inversion 

We  now  attempt  to  compare  the  upper  norm  bounds  in 
(4.9)  and  (4.39).  Since  IIa"1^  <  ||  A~ 1 1|  ^+1  E^H  ^  ,  and 
ll<5~1|l00  1  II  6e  1|l  oo+IE41ooj  lt  follows  from  (4.6)  and  (4.8) 
respectively,  that 

,  „  „  ^Sa(2-005  n2+n3)+||E  ||co}||&;1||Ji;1||oo 

(•4.40)  E  J  <  - - - , — J: - 2 - -2 -  , 

l-{EgA(2.005  n2+n3)+||E1||to}||A;1||ai 

and 

..  ,,x  „  „  t  (2.005  n2+n3)+||EJ}|«;1|J«;1|L 

Ei*  ”  ”  l-{egfi(2..005  n2+n3)+||E3||oo>||  <S-1||oo 

where  ||  E^H  ^  and  ||  E  ^  ||  ^  are  defined  by  (4.5)  and  (4.7), 
respectively.  If  computed  inverse  A  ^  and  6  ^  of  A 
and  6,  respectively  are  reasonable  approximation  to  the 
true  inverse,  we  may  assume  that 

(4.42)  {EgA(2.005  n2+n3)  +  |E1||co}||A;1|co  <  0.1  , 
and 

(4.43)  {eg6(2.005  n2+n3 )+|| E3||  ^>11  6^||  ^  <  0.1  . 

It  follows  from  (4.40)  and  (4.4l)  that 

(4.44)  ||E2||oo  <  1 . 12 { egA (  2.005  n2+n3 )+|| E1#  J\\  A;1!  J  A^fl „  , 


' 


' 


. 


. 
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and 

(4.45)  ||E4|| ^  <  1.12{Eg6(2.005  n2+n3)+||E3||oo}||6;1|| JS^IL  . 

Without  loss  of  generality,  we  can  assume  that  elements 

r. .  of  matrix  R  have  been  scaled  such  that 
ij 

(4.46)  max  I r . . I  <  1  . 

1  ij  1  - 

Therefore 

(4.47)  II  A||  „  <  n  , 

and 

(4.48)  |l B|| w  <  n  . 

Using  (4.47),  (4.48)  in  (4.5)  and  (4.7),  respectively, 

(4.49)  || E ill  m  <  e[n+n2{gA(2.005  n2+n3 ) ||  A”1!  a+  ( 2n+l)+n(n+2 ) £ } 

-  iia;4l]  . 
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and 

(*••50)  ||E3||oo  <  e  [n+n2  {gg(  2 .005  n2+n3)  j|  B"1!!  <x>+(  2n+l) 

+  n(n+2)£}|!B;1||oo]  . 


Therefore , 


(4.51)  ||E1||oo+||E3||oo  <  e[2n+n2{(2„005  n2+n3) 

x  (sAll  Ae1H  coll  A”1!!  00+gBll  B"1!!  Jl  B^IL) 

+  (2n+l)  (||Aci]|oo+||Bcxj|oo)}]  , 

2 

with  terms  of  order  e  ignored.  In  order  to  guarantee 
that  the  computed  inverse  of  A  and  6  is  a  reasonable 
approximation  to  the  true  inverse 3  we  assume  that 
II Ol-  i  1-01llAe1lh  and  ll«o1L  5  1-01||«"1|L  •  since 
l|A”1|L  i  |lRe1|L  and  !l  6  e 1  ‘I  ~  i  !lRell[h  ’  we  have  from 

(4.44),  (4.45)  and  (4.51) 


(4.9) , 


that 


’oo.: 


, 


).  1  »  Xa  '  ■'  -  o"i:  %  I 


6  6 


( 5 ) 

(4.52)  l|R"1-R”1|L  <  1.12e[(ga+g6)(2.005  n2+n3)+2n 

+  1 . 01{n2  (  2 . 005  n2+n3)(gA||A"1||£<>+gB!|B"1||co) 

+  (  2n+l )  ( |t  A"1!)  ,^+H  B"  1||  ^ ]  ||  R"1!!  J  R"1!!  „  , 

-1  -1  (S) 

where  j|  R  -R  j|  denotes  the  °°-norm  of  the  matrix 

II  Q  0  II  00 

of  round-off  errors  incurred  in  computing  the  inverse 

of  a  block-symmetric  matrix  R  using  Schur's  Algorithm  III 

( see  Section  4.1). 

Similarly,  assuming  that 

C4.53)  e(l|A|L+|B|ioo)||p;1||oo  <  0.1  , 

and 

(4.54)  2e  ( ||  A||  w+2||  Bj|  j|  Q~  1ji  w  <  0.1  , 

it  follows  from  (4.31),  (4.36),  (4.39),  (4.47)  and 
(4.48)  that 


' 


■ 


I 
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(4.55)  iir^-r;1,^ 


(C) 


;  e(||P^1||00+ll  Q“1|L)+1 . 12e 


x  {(2.005  n2+n3)(gp||p;1||  J|  P^IL+SSqII  QjfiL 


x  llQe1IL)+2n(||Pe1|L+||Qe1||  JJ  , 


-1  -1  (C) 

where  ||  R  -R  II  denotes  the  “-norm  of  the  matrix 

m  c  e  11 00 

of  round-off  errors  incurred  in  computing  the  inverse  of  a 
block-symmetric  matrix  R  using  Charmonman's  Algorithm  I 
(see  Section  4.3). 


II  Q 


Assuming  that  ||PC1||00  <  1 . 0l||  Pe1||  ^  and 

;1IL  <  1  •  01ll  Qg^ll  °°  9  we  have  from  ^*55), 


(4.56)  | R  1-R  1„ 

11  c  e  " 00 


(C) 


<  1.01  e(||Pe1||co+lQe1||o0)+1.12e 


2  2 

X  {1.01(2.005  n2+n3)(gp||P;1|i0ot2gQ|Q-1||co  ) 


+  2n(|P;1|L+l|Qe1«!)} 


. 


' 
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Since  R  ^  = 

E  i  F  ~ 
1 

_ l__  _ 

1 

,  it  follows  from  the  definition  of 

F  I  E 

“-norm  that 

(4.57) 

II  EL  i  IIR^IL  . 

and 

(4.58) 

I|F|L  i  II  Re1IL  • 

-1 

-1 

Further,  since 

P  = 

E+F  and  Q  =  E-F,  it  follows  from 

e 

e 

(4.57)  and  (4.5&) 

that 

(4.59) 

llfflL  i  2 ll  Re1iL  . 

and 

(4.60) 

IIQ^IL  i  2II Re1IL  • 

Hence,  from  (4.55), 

(4.59)  and  (4.60)  , 

(4.61)  |Rc1-r; 

,  (C) 

;  lie 

<  £[1.12(4.04(2.005  n2+n3) (gp+2gQ) 

+  32n}  +  2.02]||Re1||co||Re1|L  , 


* 


yd* 
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provided 


(4.62) 


>  1  . 


It  is  observed  in  practice  that  if  the  elements  of  a  non¬ 
singular  matrix  are  bounded  by  unity,  then  the  elements 
of  its  inverse  are  generally  greater  than  unity.  Therefore 
in  the  light  of  (4.46),  the  assumption  in  (4.62)  is 
generally  true.  It  is  thus  obvious  by  comparison  of 
(4.6l)  to  (4.52)  that  the  norm  bound  for  the  matrix  of 
round-off  errors  incurred  in  computing  the  inverse  of  a 
block-symmetric  matrix  using  Charmonman's  Algorithm  I 
(see  Section  4.3)  is  considerably  smaller  than  that 
using  Schur's  Algorithm  III  (see  Section  4.1).  In 
particular,  let  g  =  max(gA  ,gfi  ,  gQ ,  ||  A"^  oo,  ||  B"1!  J  ,  then  for 
large  _n,  the  ratio  of  the  upper  norm  bounds  in  (4,52) 
and  (4.6l)  is 


(4.63) 


Although,  the  actual  errors  depend  on  the  nature  of  the 
elements  of  the  matrices,  it  is  reasonable  and  standard 


. 


' 


. 


8t  (Sc  .t1)  ni  noliq.nu  as  »rii  «?*’•'  '~X  J  II 
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. 


70 


practice  to  say  that,  in  general,  a  method  with  smaller 
error  bound  is  better  than  that  with  larger  error  bound. 

4 . 6  Solution  of  System  of  Equations 

The  solution  of  the  system  of  linear  algebraic 
equations 

(4.64)  Rx  =  b  , 

where  the  coefficient  matrix  R  is  block-symmetric, 
may  be  computed  according  to  the  algorithm  given  by 
Charmonman  [1].  This  algorithm  requires  smaller  computer 
storage  and  smaller  number  of  multiplications.  It  is 
based  on  the  following  theorem: 

Theorem  4.5  The  solution  of  system  of  equations 
(4.64}  which  can  be  written  in  the  form 


'  A 

B 

j 

| 

J 

A 

_ 

fz_ 

is  x i  =  0.5(y^+y^)  and  x %  =  0.5(y^-y^)  where  y ^ 
and  y ^  are  solutions  of  ( A+B ) y ^  =  b  ^+b  ^ 

(A-B)y^  =  b y-b respectively . 


and 


- 


. 

. 


■ 
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The  theorem  can  be  easily  proved,  by  equating  terms  on 
both  sides  of  (4.65). 

Charmonman ' s  Algorithm  II  To  compute  the 
solution  of  (4.64): 


Step 

1 : 

Compute 

P=A+B  . 

Step 

2  : 

Compute 

Q=P-2B  . 

Step 

3: 

Compute 

Cl=bl+b2  * 

Step 

4: 

Compute 

c2=c1-2b2  . 

Step 

5: 

Compute 

y1  from  Py1=c1 

Step 

6: 

Compute 

y  2  from  Qy2=c2 

Step 

7: 

Compute 

x1=0 . 5(yx+y2)  . 

Step 

8: 

Compute 

x2=xl“y2  * 

4 . 7  An  Error  Analysis  of  Solution  of  Equations 

We  now  proceed  to  place  norm  bounds  on  the  round-off 
error  incurred  in  the  computation  of  the  solution  of 
(4.64)  using  Charmonman' s  Algorithm  II  (see  Section  4.6). 
We  will  use  the  same  notations  as  in  Section  3-2. 

4.7.1  Bounds  on  HeJ^,  ||E2|oo,  ||E3||oo  and  ||  E  4 1|  ^ 

Upper  norm  bounds  on  the  round-off  errors  incurred  in  the 
first  four-steps  of  Charmonman ' s  Algorithm  II  follow 
from  -  (2.27)  .  Thus 


■ 


' 
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(4.66) 

!lEllL  i  e(||A|L+||B||J  , 

(4.67) 

IIe2II»  5  2e(||A||oo+2||B||oo)  , 

(4.68) 

IIe3!I.  i  sdhxIL+llb^L)  , 

and 


(4.69) 

llE4ioo  ;  2e(  !  bj  „+2||  b2||J  . 

OJ 

• 

t— 

• 

Bounds  on  ||E5||oo  and  1 E g |) ^  Let  y^15  and 

denote  the  exact  and  computed  solution  of  Py=c-,, 

^  _L  JL 

then  it  can  be  easily  shown  that  the  computed  solution 
is  an  exact  solution  of  the  perturbed  system 


(4.70) 

(P+E1+K)y^1)  =  Ce1J+E3  ’ 

where  c^1'*  denotes  the  exact  value  of  cn  „  A  bound  on 


K  follows 

from  (3.68) . 

We  have 

(4.71) 

IlKlL  <  egp(2.005  n2+n3)  . 

(4.71) 


noi^JJloa  bscrjjqi  oo  9<i3  - nworft  ^  >*?  9  1  ^  I 

■ 
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Let  h  =  y^-y^,  then  from  (4.70), 


(4.72) 

|h|  =  |p;1{E3-(E1+K)y^,1)}|  , 

or 


(4.73) 

|yc1)'ye1,l  =  lpe1{E3-(El+K)yc1)}l  > 

or 


(4.74) 

l|yc1)-ye1)ll-  i  llPe1|Lt|E3IL+(||E1|co+||K|oo)||y'1)||oo} 

Using 

(4.66),  (4.68)  and  (4.71)  in  (4.74)  and 

noting 

that  E^y ^ ^ -y ^ ^ ^ ,  we  derive 

3)  0  0 

(4.75) 

IIe5IL  i  e  C  (gp (2.005  n2+n3)+|| A||  ^+1  B||  ^>11  y^1 4  „ 

+  Hb1IL+llb2LJIIpeiIL  > 

where 

gp  is  the  element  of  maximum  modulus  at  all  stages 

in  the 

reduction  of  computed  P  by  use  of  Gaussian 

elimination . 


v 


. 

' 


- 

. 


Similarly,  the  upper  norm  bound  for  the  matrix  of  round¬ 
off  errors  incurred  in  computing  the  solution  of  the 
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system  Qy^c^  is  given  by 

0.76)  II  Eg||  K  <  e  [  (gQ(  2 . 005  n2+n3)+2(||A||co+2||B|U)}||y^2)||oo 

+  2(«»>1«.+2|b2iL)]|Q;1iL  . 

where  gn  denotes  the  element  of  maximum  magnitude  in 

the  transformation  of  computed  Q  using  Gaussian  elimination 

(  2 ) 

and  y^  denotes  the  computed  solution  of  ^2~Q'2  * 

4.7.3  Bounds  on  jjE^lj^  and  j|  Eg(|  The  computational 
equation  for  the  calculation  of  x^  is 

(4.77)  x1  i  0.5{y,i1J+E5+y^2)+E6}+E71  . 

Therefore 

(4.78)  |E?|  =  |0.5(E5+Eg)+E71|  , 


or 


■ 


' 


' 
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(4.79) 

|e7IL  <  o.5(||e5||j-||e6||j+||e71il  . 

Also , 


(4.80) 

llE7ilL  ;  °-5(||E5|| „+||e6||oo)  . 

Hence , 

it  follows  from  (4.75),  (4.76),  (4.79)  and 

(4. 80) 

that 

(4.81) 

I|e7IL  <  £[(llb1IL+«b2iL)||p;1iL+2(||b:L||oo+2||b2iL)||Q;:LiL 

+  (gp(  2 . 005  n2+n3)+||A||oo+[|B||oo>||p;1||J|y'1)|.oo 

+  {gQ(2.005  n2+n3)  +  2(||A||oo+2i|B||co)}!|Q;1||  . 

Similarly,  we  have 


(4.82) 

I|s8IL  i  IIe6IL+IIe7IL+1|e8i|L  . 

where 

(4.83) 

Eg1  is  bounded  by 

l|E8lIL  <  E||xJ,1;-y'2)|L  • 

(4.83) 


r<  .'J^-Klljaipe.o  >  jlxn ' 


■ 


. 
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The  result  of  this  section  is  summarized  in  the  following 
theorem : 


Theorem  4,6  If  the  solution  of  (4.64)  is 
computed  using  Charmonman’s  Algorithm  II  (see  Section  4.6)  3 
then  the  round-off  error  in  the  computed  solution  x 
is  hounded  by 

(4.84)  ||®0-*,il  _<  - 

where  E Q  and  E ?  are  defined  by  (4.76)  and  (4.81)  3 
re spe ctively 3  and  denotes  the  exact  solution  of 

(4.64) . 

4 . 8  Comparison  of  Upper  Norm  Bounds  in  System  of  Equations 

We  now  proceed  to  compare  the  upper  norm  bounds  on 
the  round-off  errors  incurred  in  computing  the  solution  of 
system  of  equations  (4.64)  using  Charmonman’s  Algorithm  II 
(see  Section  4.6)  and  that  obtained  by  using  Schur’s 
Algorithm  II  (see  Section  3.5)  for  a  block-symmetric  co¬ 
efficient  matrix.  Norm  bound  for  the  error  in  the  computed 
x1  and  x^  can  be  obtained  from  (3.66),  (3.67),  (3-70), 
(3*72),  (3.84),  (3.85)  and  (3.86)  by  replacing  r  by  n, 

C  by  B  and  D  by  A. 


' 


' 

* 


' 
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Thus 

(4.35)  llxc1)_xe1)U<»  -  e[gB(2.005  n2+n3 )||  x^1  J 1| ^+11  b2||  M 

+  (n+l+ne  )||x^  2)||  J|  A||  J||  B^U  _+||  x< 2  > -x^2)||  J  A||  J  B^U  „  , 

and 

(1).86)  ||x^2)-x^2)||oo  <  e[||b2i|oo+||A||oo||x^2)||TO+{gA(2.005  n2+n3||  A^l 

+  l+2n+n(  n+2 )  £ }  (||  bj  m+\\  Bj|  J  x^2 )  ||  J  ||  B||  J  A^H  m 
+  gA(2.005  n2  +  n3)||x'2)||co]||A;1|L  , 

-1  /-Is 

where  A  denotes  the  exact  inverse  of  (A-BA  B). 
e 

If  we  assume  that 

(4.87)  l|B_1||  J|A||„  >  1  , 

then  the  norm  bound  for  the  round-off  error  in  computed 

x  is  given  by 
c 


' 


(a  "AS-A)  'to  38T9Vni  J0BX3  9rf3  Ef  s‘-  019riw^ 
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(4.88) 

ilx0-xeL  1  e[gB(2.005  n2+n3)||x^1)j|oo+||b2!|co 

where 

+  (n+l+n£)||x^2)j| J|a||co]!|b"1||oo+||x^2)-x^2)!|co  , 

||xf  )-f  2)|l  is  bounded  by  (4.86). 

Without  loss  of  generality,  we  can  assume  that  elements 
of  the  coefficient  matrix  R  have  been  scaled  such  that 

(4.89)  | r „ , |  <  1  . 


With  this 

assumption,  we  certainly  have 

(4.90) 

> 

8 

1  A 

V* 

(4.91) 

ll  B  |  <  n  . 

11  OO  _ 

Substituting  (4.86)  in  (4.88)  and  using  (4.90)  and 
(4.91),  it  follows  that  the  norm  bound  on  round-off  errors 
incurred  in  computing  the  solution  of  (4,64)  by  Schur's 
Algorithm  II  (see  Section  3.5)  is  given  by 


' 
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(4.92)  !|xo-xe||iS)  <  e  [gB  (  2 . 005  n2+n3)||x£1,||eo+||b2||oi> 

+  n(n+l+n£)||x^2)||a)]||B;1||(!o+En[||b2l|oo+n||x^2)||co 
+  n{gA(2.005  n2+n3)||A”1||co+l+2n+n(n+2)£} 

x  (||b1IL+n||x'2)||co)||A;1||<o+gA(2.00  5  n2+n3)||  x<2)||J 

x  ||b;1||J'a;1il  • 

Similarly,  substituting  (4.76)  and  (4.8l)  in  (4.84) 
and  using  (4.90)  and  (4.91),  we  derive  the  following 
norm  bound  for  the  round-off  errors  incurred  in  computing 
the  solution  of  (4.64)  using  Charmonman’s  Algorithm  II 
(see  Section  4.6): 

(4.93)  l|x0-xel|^C)  <  e[||x^1)||co+||y^2)j|otj+(||b1]|0(>+||b2||m) 

x  l|P;1L+2(l|b1IL+2||b2||co)||Q;1||0o 

+  ( 2n+gp( 2.005  n2+n3)}||p;1||J|y^1,||co 

+  3(6n+gQ(2.005  n2+n3)}||Q"1||co||y^2)!|co]  . 


* 
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It  is  thus  obvious  from  (4.92)  and  (4.93)  that  the 
dominating  term  in  the  expression  for  round-off  error 
incurred  in  computing  the  solution  of  (4.64)  using 
Schur's  Algorithm  II  (see  Section  3*5)  is  O(n^),  whereas 
the  corresponding  error  term  using  Charmonman’s  Algorithm 
II  (see  Section  4.6)  is  O(n^).  Wilkinson  [17]  has 
shown  that  the  round-off  error  incurred  in  computing  the 
solution  of  (4.64)  using  Gaussian  elimination  only,  is 
also  of  order  n  .  However,  the  actual  errors  depend 
on  the  nature  of  the  elements  of  the  coefficient  matrices 
considered. 


' 


' 


. 
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CHAPTER  V 


BAND  MATRICES  USING  MINIMUM  STORAGE 

5 . 1  Solution  of  Equations 

Gaussian  elimination  method  (see  Section  2.2)  can 
be  used  to  solve  a  system  of  linear  equations  with  large 
band  coefficient  matrix.  We  will  consider  a  general  band 
matrix  in  which  the  number  of  co-diagonals  on  each  side  of 
the  main  diagonal  are  not  equal.  The  concept  of  band 
matrix  is  useful  only  when  the  band  width  is  appreciably 
smaller  than  the  order  of  the  matrix.  The  band  width 
depends  upon  the  ordering  of  the  unknowns  and  the  equations 
in  the  system.  Evidently,  row  or  column  interchange  will 
result  in  increase  of  band  width  and  accordingly  more 
storage  space  is  required.  Further,  if  complete  pivoting 
is  used  the  band  structure  may  be  completely  lost. 

However,  with  partial  pivoting,  the  number  of  upper  co¬ 
diagonals  may  increase  but  the  maximum  increase  can  not  be 
more  than  the  number  of  lower  co-diagonals.  Also,  the 
band  characteristics  of  the  given  matrix  are  preserved. 

In  order  to  take  advantage  of  the  band  form  of  the  given 
coefficient  matrix,  only  the  elements  on  the  band  are 
stored  in  the  computer.  One  way  to  accomplish  this  is 


. 

• 

' 

. 
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to  store  the  elements  on  the  band  row-wise.  The  LU 
decomposition  of  the  given  band  matrix  is  then  performed 
using  an  algorithm  given  in  [7]  based  on  Gaussian  elimina¬ 
tion,  which  uses  minimum  storage.  Since  partial  pivoting 
is  used  the  upper-triangular  matrix  U  will  involve  as 
much  storage  as  A  in  the  worst  case.  Therefore, 
additional  storage  must  be  provided  for  the  lower-triangular 
matrix  L  or  if  L  is  not  saved,  the  right-hand  sides 
must  be  processed  simultaneously. 

The  algebraic  basis  of  the  decomposition  of  band 
matrix  A  into  lower  and  upper  triangular  matrices  is 
contained  in  the  following  theorem: 


Theorem  5.1  If  cl  square  band  matrix  of  order 


n  with  p  lower  and  q  upper  oo-diagonals  has  an  LU 

decomposition ,  then  L= ( £  .  °)  and  U= ( u .  .)  are  triangular 

I'O  'c  C 

band  matrices .  That  is , 


7*  0  only  for  i=j ,  , 


•  *  •  3 


C+V  , 


z+q  . 


The  proof  of  this  theorem  follows  from  Theorem  2.1. 


' 
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5 . 2  Upper  Bound  for  the  Pivot  in  the  rth  Reduced 

Band  Matrix 

Let  A  denote  an  unsymmetric  non-singular  square 

band  matrix  of  order  n  with  p  lower  and  q  upper 

co-diagonals.  Then  using  Gaussian  forward  elimination, 

the  original  matrix  A=A^^  is  transformed  successively 

into  matrices  A^2\  A^  \  .  .  .  ,A^r\ .  .  .  sA^n^  ,  such  that 
( r )  ( 1 ) 

each  A  is  equivalent  to  A  and  the  final  reduced 

matrix  A^n^  is  upper-triangular.  We  now  proceed  to  place 

an  upper  bound  on  the  magnitude  of  pivotal  element  in  the 
t  h  ( r  ) 

r  reduced  matrix  A  using  partial  pivoting.  In  the 

first  few  stages  of  reduction  of  A,  the  number  of  elements 
in  the  pivotal  column  Is  p+1  and  the  corresponding  rows 
to  be  transformed  are  p  in  number  excluding  the  pivotal 
row.  These  rows  include  one  row  which  has  not  yet  been 
modified.  Thus  from  the  growth  of  the  elements  in  the  pivotal 
sequence,  it  follows  that  every  successive  pivotal  element, 
in  the  worst  case,  is  the  sum  of  preceding  p<n-l  elements 
plus  a,  where  |a|^|<a.  The  undefined  elements  in  the 
pivotal  sequence  are  treated  as  zero. 


Therefore , 


<  t 

-  1  r-p 


+  1 1  (  x  +.  .  .  +  t  +a 

1  r-(p-l)  1  1  r-1 1 


(5.1) 


3 


■ 

' 


■ 
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where  tr(r=p+l ,p+2  , .  .  .  ,n)  denote  the  pivot  at  the  r ^ 
transformation.  Since  each  t  is  bounded  by  2r_1a, 
we  have  from  (5.1)  s 

(5.2)  jt^|  <  {2r  1(l-2  p)+l}a,  r=p+l 5p+2 . . . ,n, 

p<n-l  . 


The  following  theorem  is  a  result  of  (5.2): 


Theorem  5 . 2  Given  a  non- singular  band  matrix 
A  of  order  n  with  p  lower  and  q  upper  eo -diag onals 3 


then  using  Gaussian  elimination  with  partial  pivoting 3 

r 

, r-1  „  „ 

.  .  .  j  p  3 


(5.3) 


a 


(r ) 

'id 


r-1 

2  a j  r=l  3 2  3 


<  \ 


{2r  1 ( 1-2  V)+l}a3  r=p+l 3p+2. . . 3n3 

p<n- 1 3 


where 


(5.4) 


(1)  s  , 

a  ,  o  <  a  , 

1  - 


and  denotes  the  (i*j)  element  of  the  k  ^  reduced 

matrix  A  ^  with  A  ^’=A. 


' 

' 


' 
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5 . 3  An  Error  Analysis 

In  this  section  we  place  a  norm  bound  on  the  round¬ 
off  error  incurred  in  the  solution  of  a  system  of  linear 
algebraic  equations 

(5.5)  Ax  =  b  } 

with  non-singular  band  coefficient  matrix  of  order  n 
using  Gaussian  elimination  with  partial  pivoting.  We 
assume  that  A  is  unsymmetr'ic  with  p  lower  and  q 
upper  co-diagonals.  We  may  write  the  coefficient  matrix 
as 
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It  has  been  shown  (for  example  see  Pox  [5])  that  the  process 

of  forward  elimination  with  partial  pivoting  is  equivalent 

to  successive  premultiplication  of  (5.6)  by  matrices 

L-,  PL  ,L~R0  ,  .  .  .  ,L  ,  R  ,  where  L.  are  unit  lower-triangular 
1  1J  d  d  n-1  n-1  i  ° 

and  R,  denotes  the  permutation  matrix  R. .  which  accom- 
1  ij 

plishes  the  partial  pivoting  at  the  ith  stage.  Since 
the  interchange  of  rows  does  not  involve  round-off  error, 
we  can  assume,  without  loss  of  generality  that  the  co¬ 
efficient  matrix  has  been  reordered  such  that  no  interchange 
of  rows  is  necessary.  Then  it  can  be  easily  shown  that 
each  L„  is  of  the  form 

l 


0  1 
0  0  1 


(5.7)  Lj_ 


0  0  0 
0  0  0 
0  0  0 


0  0  0  0  •  °  0  mp+i,i 

0  0  o  . .  •  0  0  0 


1 

0 


1 


0  m„ .  1 
1+1,1 


1 

0 


' 


■ 


■ 


■  1  2  "iL:' 


‘ 

■ 


. 
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for  i  <  n-p ,  and 


i 


(5.8)  L±  = 


0 


0 


0 


0 

0 


0  . 


0  •  o  o 

0 


1 


0  ...  0 


1 


0  ...  0  m. .  1 

1+1,1 


0  m  .  0 

n-1 ,  i 


0  m 


n,i 


0 


1 

0  1 


for  n-p  <  i  <  n-1,  where 


(5.9) 


a 


(k) 

Ik 


IkT  5 

kk 


k  =  1 , . . . ,  n-p , 

1  =  k+1 , . . . ,p+k , 


and 


(5.10) 


m0 , 
lk 


(k) 

aik 

kk 


k  =  n-p+1 , . . . ,n-l , 
i  =  k+1 , . . . ,n . 


. 


. 


. 


' 


■ 
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Further,  it  can  be  shown  that  matrices  J=L  ,L  0...Ln 
5  n-1  n-2  i 

and  .  .  .  Ln^  are  unit  lower-triangular,  and 

J  and  L  have  the  form 


1 


m 


21 


1 


m 


31 


m 


32 


(5.11)  J  = 


m  ,  n  _  m  ,,  n 

p+1,1  P  +  1 j  2 


0 

0 


m 


P+2,2 

0 


0 


0 


1 


m 


P+1,3 


m 


P+2,3 


m 


P+3,3 


0  m  o  «  ©  m  1 

n,n-p  n,n-l 


and 


, 


. 


' 
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(5.12)  L  = 


1 


-m 


21 


1 


-m^  -m32  1 


'mp+l,l  mp  + 1  a  2  mp  +  l , 3 


0 

0 


-m 


P+2,2 
0 


0 


0 


-m 


P+2,3 


e  o  e 


-m 


P+3,3 


0  -m 


n,n-p 


-m 


n  ,n-l 


We  have  in  fact  performed  a  process  equivalent  to  the 

decomposition  of  A  into  lower  and  upper  triangular 

matrices,  the  upper-triangular  matrix  being  the  final 

( n ) 

reduced  matrix  A  which  may  have  up  to  p+q  co-diagonals. 

We  may  therefore  write 


(5.13) 


A 


or 


9riJ  ,7:  v  '  1  S  'I9CI  ''  91 
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(5.14) 


A  =  LU  , 


in  which  L,  the  lower  triangular  matrix  is  given  by  (5.12). 
The  equality  sign  in  (5.14)  can  not  be  expected  to  be 
exactly  satisfied  due  to  round-off  errors  incurred  in 
various  computations  performed.  Wilkinson  [16],  has  shown 
that  computed  L  and  U  satisfy  the  perturbed  equation 

(5.15)  LU  e  A+E  , 

where  E  is  the  perturbation  matrix.  The  solution  of 
(5.5)  is  obtained  by  solving  the  triangular  systems 

(5.16)  Ly  =  b  , 


and 

(5.17)  Ux  =  y  . 

The  computational  equations  for  computed  y  and  x  in 
(5.16)  and  (5-17)  are  given  by 

(5.18)  (L+6L )y  e  b  , 


and 


' 


. 


' 


' 
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(5.19)  (U+6U)x  E  y  . 

Hence,  computed  x  satisfies  the  computational  equation 

(5.20)  (L+6L) ( U+6U) x  =  b  . 

From  (5.15), 

(5.21)  (A+E+L6U+6LU+5L6U)x  E  b  . 

If  x  and  x  denote  the  exact  and  computed  solution  of 
0  0 

(5.5),  then  since  Axg=b 3  it  follows  from  (5.21)  that 

(5.22)  ||x0-xe|U  <  II  A-1!!  „  C II  E||  „+l|  L||  J|  6U||  „+[[  6L||  J|  U||  „ 

+  ll«Lll  J5uIL)llxcllc»  • 

5.3.1  A  Bound  for  IIeII^  Using  technqiues  similar  to 

those  employed  by  Wilkinson  [19  ]it  can  be  easily  shown  that, 

in  general  for  an  unsymmetric  band  matrix  A  with  p 

lower  and  q  upper  co-diagonals,  the  perturbation  matrix 

E  =  (e. .)  is  defined  by, 

J 


■  c 


■ 


* 


_ 

' 


. 


. 


I 


- 
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i=2,...,p+l,  and  j=2,...,q+l  ; 
i=2,...,p+l,  and  j =q+2  , . . . ,q+i  ; 
i=2j...,p+l,  and  ; 

i=p+2  , . . . ,n-p-q ,  j  =i-p  , . . . ,p+q+i  ; 
i=n-p-q+l, . . . ,n,  and  j=i-p,...,n; 

where  the  element  re , ( r=l , 2  , . .  .  ,p )  represents  the  number 
of  times  the  corresponding  element  in  the  matrix  A  was 
modified  during  the  process  of  reduction  of  matrix  A  into 
upper-triangular  matrix  using  Gaussian  elimination, 

andN 

l 

(5-24)  |  r’e  |  <  3gre  ,  j  >  i  , 

and 

(5.25)  | re |  <  g(3r-2)e  ,  i  >  j  . 


(i-)e , 
(i-j+q+l)e, 
j  e  3 

/ 

ep+i,j+p+i-i’ 

ep+i,j+p+i-i5 

0  3  otherwise , 
V 


g  is  defined  as 


. 


. 

* 

* 

bne 


9 4 


(5.26) 


max  |  (k) 

aij 


A  bound  on  g  has  already  been  given  in  Section  5.2.  It 
can  be  easily  seen  that  starting  from  the  second  row 
through  (p+l)^*1  row,  the  number  of  non-zero  elements 
in  the  i^*1  row  of  the  perturbation  matrix  E  defined  by 

(5.23)  is  p+q+i-1,  (i=2 , . . .  ,p+l)  ;  thereafter  in  each 
successive  row  this  number  decreases  by  1.  Thus  2p+q  is 
the  maximum  number  of  non-zero  elements  in  any  row  of 

(5.23) .  The  result  of  this  section  is  summarized  in  the 
following  theorem: 


Theorem  5.3  Given  a  non-singular  band  matrix 
A  of  order  n  with  p  lower  and  q  upper  oo-diag onals . 
Then  if  2p+q<n ,  the  matrix  E  of  round-off  errors  incurred 
in  decomposing  A  into  lower  and  upper  triangular  matrices 
using  Gaussian  elimination  with  partial  pivoting  is  bounded 
by 

(5.27)  ||ff|L  <  zg{p<Sp+3q+l)}, 


where  g  is  defined  by  (5.26). 


. 

• 

- 


' 


. 


- 


/ 


' 


. 
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5.3.2  Bounds  on  j|  L||  ^  and  ||  U||  ^  Since  partial 
pivoting  ensures  that  the  multipliers  lmij I  f  therefore 
it  immediately  follows  from  (5.12)  that 


(5  -  28)  || L|| ^  <  (p+1)  . 

We  now  proceed  to  place  an  upper  norm  bound  on  the  upper- 
triangular  matrix  U.  The  matrix  U  is  of  the  form 


(5.29)  U 


a 


(1) 

12 


a 


(2) 

22 


1  jP+q+1 

(2)  (2) 

a2,p+q+l  a2,p+q+2 


( n-1 ) 
an-l ,n-l 


a 


(n-1) 
n-1  ,n 

(n) 

"nn 


with  p+q  upper  co-diagonals,  in  the  worst  case.  Since 

i  ( r* )  i  \ 

g  is  the  maximum  element  of  any  |A  |,  ( r=l , 2 , . . . ,n ) , 
therefore  from  (5.29)  3 


* 


* 


' 


' 

' 


, 


. 


■ 
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(5.30)  ||U|L  <  g(p+q+l)  . 

5.3.3  Bounds  for  ||  6L||  ^  and  ||  6U||  ^  Error  analysis 
of  the  solution  of  triangular  sets  of  equations  (5*16) 
and  (5.17)  is  analogous  to  that  given  by  Wilkinson  [17]. 
Using  a  technique  similar  to  that  employed  by  Wilkinson 
[17],  it  follows  using  (2.18),  (2.19)  and  (2.20)  that 
computed  solution  of  (5.16)  is  an  exact  solution  of 
(L+6L)x=b,  where  6L 


is  bounded  by 


8 


, 

' 


* 


' 


, 

' 
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Since  the  elements  rrn  .  of  lower-triangular  matrix  L 
are  bounded  by  unity,  therefore 


(5.32) 


I  5LII  oo  5  ^e(p+l)  (p+2)  . 


Similarly,  it  can  be  shown  that  the  computed  solution  of 
(5.17)  is  an  exact  solution  of  (U+6U)x=y,  where  6U 
is  bounded  by 


. 


x.  •  .  fl£»  t"- 


, 


' 


(  c  — 1 )  I  u 
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U) 

V  I 

‘O 

on 

on 

• 

Ln 
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where  c  =  p+q,  d  =  n-c  , 


and 


Hence , 

(5. 3*0  ||  6U||  ^  <  e^(p+q+l)  (p+q+2)  . 

The  following  theorem  is  a  result  of  (5.22),  (5.27), 
(5.28),  (5.30),  (5.32)  and  (5-34): 

Theorem  5.4  Let  x  and 

-  Q 

and  computed  solution  of  the  system 
equations , 

(5.35)  Ax  =  b 

where  A  is  a  band  ooeffioient  matrix  with  p  lower 

and  q  upper  co-diagonals .  Then ,  if  the  terms  of  order 
2 

e  are  neglected ,  the  round-off  error  in  computing  the 
solution  of  (5.35)  using  Gaussian  elimination  with 


x  denote  the  exact 
c 

of  linear  algebraic 


- 


' 


‘ 


' 
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partial  pivoting  is  bounded  by 


(5.36)  <  keg{2p(3p+3q+l)+(p+l) (p+q+1) 


*  (2p+q+4)}\\A-e1\\Jxo 


for  2p+q  <  w  . 

5 . ^  Inversion 

Inversion  of  band  matrices  using  minimum  storage 
can  be  carried  out  by  solving  (5.5)  with  right-hand 
side  replaced  by  an  identity  matrix  of  order  n.  The 
solution  matrix  is  stored  in  the  successive  locations  of 
right-hand  side  matrix,  therefore  no  extra  storage  is 
necessary  for  the  solution  matrix. 

5 . 5  Error  Analysis  of  Inversion 

The  error  analysis  of  inversion  of  a  band  matrix  is 
analogous  to  that  given  in  the  last  section.  Using  the 
results  of  Section  5.3  it  can  be  shown  that, 

(5.37)  ||  A”1-A“1||  ^  <  Sseg{2p(3p+3q+l)  +  (p+l)(p+q+l) 


* 


; 


' 


102 


where 
of  A, 
nitude 
lower 


X  (2p+q+l))}||Ae1||J|  A01 


and  A^  denote  the  exact  and  computed  inverse 
respectively.  g  is  the  element  of  maximum  mag- 
at  all  stages  in  the  decomposition  of  A  into 
and  upper-triangular  matrices. 
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CHAPTER  VI 

NUMERICAL  RESULTS 

Numerical  experiments  were  carried  out  on  IBM  System 
360/67  computer  at  the  University  of  Alberta.  FORTRAN  IV 
source  language  was  used  for  general  and  block-symmetric 
matrices,  and  APL\360  for  band  matrices.  Several  matrices 
with  known  exact  inverses  were  used  in  testing  the  algorithms 
for  inversion  and  solution  of  linear  algebraic  equations. 

The  right-hand  sides  of  the  systems  of  equations  were  chosen 
arbitrarily.  The  exact  solution  of  rhe  system  of  equations 
Ax=b  was  obtained  by  computing  the  product  A-1b,  where  A  1 
is  the  exact  inverse  of  the  coefficient  matrix  A,  and  b 
the  vector  of  right-hand  sides.  Norm  bounds  for  the  round¬ 
off  errors  as  given  in  Chapters  III,  IV  and  V  were  computed 
in  the  above  computer  programs .  A  comparison  is  made  between 
the  predicted  and  actual  relative  round-off  errors. 

6 . 1  Test  Data 

The  following  test  matrices  were  used  for  inversion  and 
solution  of  system  of  equations.  The  matrices  denoted  by 
A-^ ,  A^ ,  A^ ,  Ajj ,  A^  are  due  to  Newman  and  Todd  [10],  while 
matrices  Ag  and  A^  are  due  to  Charmonman  and  Julius  [2]. 


' 

' 


■ 

■ 

- 


' 
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The  matrix  A1  is  defined  by 


(6.1) 


a.  . 
ij 


-il/2 


n+1 


sin 


IjTT 

n+1 


Since  A^  is  symmetric  and  orthogonal,  we  have 


(6.2) 


A  ^  =  A 
A1  A1  ’ 


The  matrix  is  defined  by 


(6.3) 


i/j  j 

i  <  J  ; 

a.  .  =  < 

1J  1  j/i  . 

i  >  j  . 

The  inverse  of  A^  is  symmetric,  triple-diagonal 
matrix  with 


(6.4) 


a.  .  =  { 

ij 


4i- 


4i2-i  * 


n 


2n-l 

-1(1+1) 

21+1 


0 


i= j ,  i  <  n  ; 


i  =  j=n 


,  i  <  j ,  |  i-j I =1  ; 


aj  ±  ,  i  >  j  ,  I  i-j  |  =1  ; 


i-j  I  >  1  • 


. 


v 


. 


. 


' 


c  . 


•  t 


:  1  - 
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The 

matrix 

A3  is 

defined  by 

(6.5) 

aij 

=  2  min(i , j )-l 

and  its 

inverse  is  a  triple-diagonal 

1.5 

,  i=j=l  ; 

0.5 

,  i=j=n  ; 

(6.6) 

a.  .  — 

ij 

1 

,  l<i=j<n  ; 

-0.5 

j  1 i-J l=i  ; 

0 

\ 

,  | i-j | >1  . 

The 

matrix 

A4  is 

defined  by 

f 

-2  , 

i=j  ; 

(6.7) 

ii 

•H 

1  , 

• 

1 — 1 

II 

1 

•H 

0  , 

i — 1 

A 

arD 

1 

•H 

Its  inverse  is  given  by 


(  -i ( n-j  +1 ) 
n+1 


(6.8) 


a.  .  =  -■ 
ij 


a .  . 

Ji 


1  <  j 

i  >  J 


v 


,  l  ....  *  '■  v  '  '  •  5 


' 


v 
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The  matrix  Ar 


is  defined  by 

r 


(6.9) 


U 


5  , 

6  , 

-<  , 
1  , 

v  °  • 


i=j=l,  i=j=n  ; 

i=j  ; 

I  i-J l  =  i  ; 
|i-j|=2  , 
otherwise  . 


It  can  be  noted  that 


(6.10) 


a,-  ='  a; 


where  A^  is  defined  by  (6.7).  Thus 


(6.11) 


Aj1  =  (A,,1)  , 


where  aT  is  defined  by  (6.8). 

Matrices  Ag  and  A^  have  been  constructed  using 
the  following  definition  due  to  Charmonman  and  Julius  [2]. 


Definition  6.1  An  r-circulant  is  a  square  matrix 

of  order  n  in  which  the  i ^  row3  i=2333...3n3  is 

•  •  •  • 

obtained  from  the  (i-1)  row  by  cyclically  shifting  each 

element  r  -places  to  the  right. 


‘ 


' 
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The  word  "row"  may  be  replaced  by  "column"  if  the 

word  "right"  is  replaced  by  "down". 

The  matrices  and  are  r-circulants  with 

first  rows  {a , a+h , . . . ,a+ (n-1 )h}  and  {a, ah , . . . ,ahn-1} , 

respectively.  The  inverse  of  Ag  is  the  s-circulant 

T 

with  the  first  column  {b-a,  b,  .  . .  ,b,  b  +  a,}  ,  where 
s  satisfies 

(6.12)  rs  =  kn+1  , 


for  some  integer  k,  and 


(6.13) 


b 


2 

> 

n2{2a+(n-l)h} 


and 


(6.14) 


a 


1 

n  h 


The  inverse  of  A ^  is  also  an  s-circulant,  where  s  is 

T 

defined  as  in  (6.12)  with  first  column  (b  ,0 , . . . ,0 ,-hb }  , 


where 


1 

a(l-hn) 


(6.15) 


b 


' 


« 
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6 . 2  Summary  of  the  Results 

6.2.1  General  Matrices  Matrices  A^,  A^ ,  and  A^  of 
Section  6.1  were  used  in  testing  the  algorithms  given  in 
Chapter  III  for  inversion  of  general  matrices  and  solution 
of  system  of  equations  with  general  coefficient  matrices. 

Table  6.1  gives  the  predicted  relative  errors  and  the 
actual  relative  errors  incurred  in  the  computation  of 
inverses  of  general  matrices  using  Schur’s  Algorithm  I 
(see  Section  3*1)  and  Gaussian  elimination  with  complete 
pivoting.  The  predicted  abs'olute  error  ||  A_^-A-^||  using 
the  above  algorithms  is  given  by  (3.60)  and  (3.17), 
respectively,  and  the  predicted  relative  error  is  defined 
as 


Since  the  exact  inverse  of  matrix  A  is  known,  actual 
relative  error  can  be  easily  computed. 

Table  6.2  gives  the  predicted  relative  errors  and  the 
actual  relative  errors  incurred  in  the  computation  of 
solution  of  systems  of  equations  using  Schur’s  Algorithm  II 
(see  Section  3-5)  and  Gaussian  elimination  with  complete 
pivoting.  The  predicted  absolute  error  j|xc-xe||oo  using 


' 

.. 

* 

V 

xfcaj^  o  -cc  5ra‘  I  si  snag  1©  evr : 

■ 


’  ’  •- 

3 1  <  no  '.■£.  -i  >  »f ( r-  •  fl1  D-  ' 
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these  algorithms  i 

(3.17)  with  A"1 

0 

predicted  relative 

(6.17) 


s  given  by  (3.77)  and/or  (3.87),  and 
replaced  by  x  ,  respectively.  The 
error  is  defined  to  be 

c-xelb 

xelh 


The  exact  solution  can  be  computed  since  the  exact  inverse 
of  the  coefficient  matrix  is  known.  As  a  result,  therefore, 
the  actual  relative  error  can  also  be  easily  computed. 


6.2.2  Block-Symmetric  Matrices  Matrices  Ag  and 
Aj  of  Section  6.1  were  used  in  testing  the  algorithms 
given  in  Chapter  IV  for  inversion  of  block-symmetric  matrices 
and  solution  of  system  of  equations  with  block-symmetric 
coefficient  matrices.  The  three  parenthesized  numbers  after 
the  matrix  Ag  and  A^  are,  respectively,  the  values  of 
a,  h,  and  r  used  to  construct  these  matrices. 

Table  6.3  gives  the  predicted  relative  errors  and 
actual  relative  errors  incurred  in  the  computation  of 
inverses  of  block-symmetric  matrices  using  Schur's  Algorithm 
III  (see  Section  4.1),  Charmonman ' s  Algorithm  I  (see  Section 
4.4)  and  Gaussian  elimination  with  complete  pivoting.  The 
predicted  absolute  error  ||A-1-A-1||  using  the  above 
algorithms  is  given  by  (4.9),  (4.39)  and  (3.17)  respect¬ 
ively,  and  the  predicted  relative  error  is  defined  by  (6.1 6). 
The  actual  relative  error  can  be  computed  since  the  exact 


inverse  is  known. 


- 


. 


an\' 

■ 

j, 

■ 

' 
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Table  6.4  gives  the  predicted  relative  errors  and  the 
actual  relative  errors  incurred  in  the  computation  of  solu¬ 
tion  of  systems  of  equations  with  block-symmetric  coeffi¬ 
cient  matrices  using  Schur’s  Algorithm  II  (see  Section  3-5), 
Charmonman’s  Algorithm  II  (see  Section  4.6),  and  Gaussian 
elimination  with  complete  pivoting.  The  predicted  absolute 
error  ||  x  — x  ||  using  these  algorithms  is  given  by  (3.87), 

Lx 

(4.84),  and  (3.17)  with  A  ^  replaced  by  x  respectively, 

c  c 

and  the  predicted  relative  error  is  given  by  (6.17). 

The  exact  solution  can  be  computed  since  the  exact  inverse 
of  the  coefficient  matrix  is  known. 

6.2.3  Band  Matrices  Matrices  A^  and  A^  of  Section 
6.1  were  used  in  testing  the  algorithm  given  in  Chapter 
V  for  inversion  of  band  matrices  and  solution  of  system  of 
equations  with  band  coefficient  matrices. 

Table  6.5  gives  the  predicted  relative  errors  and  the 
actual  relative  error  incurred  in  the  computation  of 
inverses  of  band  matrices  using  Gaussian  elimination  with 
partial  pivoting  with  and  without  minimum  storage.  The 
predicted  absolute  error  J|  A  1 — A~  ^  is  given  by  (5.37) 
and  (3-17)  respectively,  and  the  predicted  relative 
error  is  given  by  (6.16).  The  actual  relative  error  is 
computed  using  the  given  exact  inverse. 
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Table  6.6  gives  the  predicted  relative  errors  and 
actual  relative  error  incurred  in  computing  solution  of 
systems  of  equations  with  band  coefficient  matrices  using 
Gaussian  elimination  with  partial  pivoting  with  and  with¬ 
out  minimum  storage.  The  predicted  absolute  error 
II  xc-xei|  oo  is  given  by  (5.36),  and  (3.17)  with  A”1 

replaced  by  x  ,  respectively,  and  the  predicted  relative 

c 

error  is  given  by  (6.17).  The  exact  solution  and  there¬ 
fore  the  actual  relative  error  can  be  computed  since  the 
exact  inverse  of  the  coefficient  matrix  is  known. 
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GENERAL  MATRICES-RELATIVE  ERROR 
IN  INVERSION 


1 

Matrix 

- r 

i 

Order 

Schur’s  Algorithm  I 

Gaussian  elimination 

Predicted 

Actual 

Predicted 

Actual 

Ai 

5 

0.1633E+03 

0.4730E-05 

0 . 3116E-03 

0. 8022E-06 

A  ' 

2 

I 

5 

0.1158E+04 

0 .1503E-05 

0.6680E-03 

0 . 4920E-06 

A3 

5 

0.7690E+03 

0 .2176E-05 

0.1503E-02 

0 . 9970E-06 

ai 

10 

0 . 1892E+06 

0 .9018E-02 

0 . A022E-02 

0 . 7940E-05 

A2 

10 

0 .6921E+04 

0  /4385E-05 

0.1030E-01 

0 . 2217E-05 

A3 

10 

0 .4564E+06 

0.7764E-05 

0 .2175E-01 

0 . 3412E-05 

A1 

11 

0 . 3244E+06 

0 . 8103E-01 

0 .5700E-02 

0 . 1049E-04 

A2 

11 

0.6922E+04 

0.4314E-05 

0 . 1501E-01 

0.2691E-05 

A3 

11 

0 . 4564E+06 

0.1016E-04 

0.3151E-01 

0.1555E-04 

ai 

12 

0.6150E+06 

0. 3881E-01 

0.7978E-02 

0.9466E-05 

A2 

12 

0.2790E+05 

0 . 2287E-05 

0.2116E-01 

0.3332E-05 

A3 

12 

0 . 1576E+06 

0.9605E-05 

0. 4423E-01 

0.7734E-05 

ai 

13 

0.1415E+07 

0.5737 

0 . 1101E-01 

0.9613E-05 

A2 

13 

0 . 8641E+05 

0 . 6597E-05 

0.2902E-01 

0.2588E-05 

A3 

13 

0 . 4564E+06 

0.7130E-05 

0 . 6046E-01 

0 . 7515E-05 

A1 

14 

0.1999E+07 

0.2713 

0.1431E-01 

0 . 1076E-04 

A2 

14 

0 . 8641E+05 

0.6073E-05 

0. 3889E-01 

0.3457E-05 

A3 

14 

0 . 4564E+06 

0 . 8556E-05 

|  0.8078E-01 

(  0.9014E-05 

' 
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TABLE  6.1  (Continued) 


1 

i  Matrix 

Order 

Schur's  Algorithm  I 

Gaussian  elimination 

Predicted 

Actual 

Predicted 

Actual 

A1 

15 

0 . 2850E+07 

0.3168E+01 

0.1852E-01 

0.4934E-05 

A2 

15 

0.8642E+05 

0 . 6801E-05 

0 .5108E-01 

0.4577E-05 

A3 

15 

0 . 1124E+07 

0 . 7053E-05 

0.1058 

0 . 1699E-04 

16 

0.4474E+07 

0 .1418E+02 

0 .2466E-01 

0 . 1026E-0  4 

A2 

16 

0 . 2223E+06 

0.1174E-04 

0 .6594E-01 

0 . 4782E-05 

A3 

16 

0.1124E+07 

0.1035E-04 

0.1363 

0 . 2448E-04 

ai 

17 

0.6870E+07 

0 .1191E+02 

0 . 3040E-01 

0.8671E-05 

A2 

17 

0 .5020E+06 

0 .1421E-04 

0 . 8381E-01 

0.8317E-05 

A3 

17 

0 .1124E+07 

0.1131E-04 

0.1730 

0.2553E-04 

ai 

18 

O.IOO6E+O8 

0.1096E+03 

0 . 3743E-01 

0 . 8416E-05 

A2 

18' 

0.5020E+06 

0.1559E-04 

0.1051 

0.7826E-05 

A3 

18 

0.2459E+07 

0.9578E-05 

0.2164 

0. 25HE-04 

ai 

19 

0 .1428E+08 

0.4851E+03 

0 .4652E-01 

0. 8435E-05 

A2 

19 

0.1027E+07 

0.1628E-04 

0.1302 

0 . 8280E-05 

A3 

19 

0 .4910E+07 

0 . 2604E-04 

0.2676 

0.1855E-04 

ai 

20 

0 .1970E+08 

0 . 5608E+03 

0.5745E-01 

0. 1006E-04 

A2 

20 

0.1027E+07 

0 . 2249E-04 

0.1595 

0.5574E-05 

A3 

20 

0 . 4910E+07 

0.3071E-04 

!  0.3274 

0 . 4432E-04 

• 
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TABLE  6.2  (continued) 
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BLOCK-SYMMETRIC  MATRICES-RELATIVE 
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TABLE  6,3  (continued) 
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TABLE  6.4  (continued) 
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TABLE  6  .  5 


BAND  MATRICES-RELATIVE  ERROR 
IN  INVERSION 


Matrix 

Order 

Gaussian 
elimination 
with  minimum 
storage 

Gaussian 

elimination 

Actual 

Predicted 

A4 

5 

1.6786E-14 

1.0499E-13 

2.6923E-15 

A5 

5 

4. 2002E-13 

7 . 2827E-13 

7.3275E-14 

10 

5 .1292E-14 

2 . 199  2E-12 

3 . 4681E-14 

A5 

10 

5.6658E-12 

6.7345E-11 

1.9081E-12 

A4 

11 

6.1042E-14 

3.4306E-12 

5 . 0113E-14 

A5 

11 

8.2035E-12 

1.2781E-10 

3.9853E-12 

A4 

12 

7. 0721E-14 

5 . 0937E-12 

6.8709E-14 

A5 

12 

1.3292E-11 

2  6  2620E-10 

7.6360E-12 

A4 

13 

8. 2019E-14 

7.4282E-12 

6 . 2204E-14 

A5 

13 

1.5469E-11 

3. 8841E-9 

1 . 2 327E-11 

a4 

14 

9. 3259E-14 

1.0448E-11 

1 . 7761E-14 

a5 

14 

2.0234E-11 

6 . 3404E-10 

3 . 0414E-11 

a4 

15 

1.0611E-13 

1 , 4499E-11 

1 . 517  4E-14 

A5 

! 

15 

2 . 6719E-11 

1 . 0122E-9 

3.9568E-12 

A4 

16 

1.1890E-13 

1.9574E-11 

1. 8581E-14 

A5 

16 

3 . 4077E-11 

1 . 5552E-9 

4.2189E-12 
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TABLE  6.5  (continued) 


1 

\ 

! 

Matrix  j 

Order 

i 

; 

t 

Gaussian 
elimination 
with  minimum 
storage 

Gaussian 

elimination 

Actual 

Predicted 

1 

A4 

17 

1.3331E-13 

! 

2.6149E-11 

2.3157E-14 

a5 

17 

I 

4.3193E-11 

2 . 34  88E-9 

S-SBATE-ia 

a4 

18 

1.4766E-13 

3 . 4181E-11 

2.7709E-14 

A5 

18 

5 . 364  3E-11 

3. 4426E-9 

4.6286E-11 

a4 

19 

1 . 6361E-13 

4  0  4308E-11 

3.3369E-14 

A5 

19 

6 . 6299E-11 

4.9776E-9 

2.6218E-11 

a4 

20 

1.7952E-13 

5.6434E-11 

3.9568E-14 

A5 

ft 

20 

8.0606E-11 

7 . 0247E-9 

5 .40792-11 
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CHAPTER  VII 

CONCLUSIONS  AND  SUGGESTIONS 
FOR  FUTURE  RESEARCH 


7  •  1  Conclusions 

The  following  conclusions  are  drawn  from  the  numeri¬ 
cal  results : 

1.  As  can  be  expected,  the  predicted  round-off 
error  incurred  in  inverting  matrices  and  solving  systems  of 
equations  using  the  algorithms  described  in  Chapter  III  are 
substantially  larger  than  the  actual  round-off  error.  This 
is  due  to  the  fact  that  the  theoretical  results  are  obtained 
after  allowing  for  the  maximum  round-off  error  in  each  cal¬ 
culation  and  no  allowance  is  made  for  the  statistical  effect. 
The  results  obtained  in  Chapter  III  therefore  represent 
upper  bounds  which  are,  in  general,  much  higher  than  the 
actual  errors.  It  is  important,  however,  to  note  that  cer¬ 
tain  matrices  can  be  constructed  such  that  the  actual  round¬ 
off  error  Is  close  to  the  predicted  bound, 

2.  The  upper  bounds  for  the  round-off  error  incurred 
in  computing  the  inverse  and  solution  of  a  system  of 
equations  using  Schur's  algorithms,  described  in  Chapter  III, 
are  comparatively  larger  than  the  corresponding  errors  using 


1 

. 
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Gaussian  elimination.  This  supports  the  remark  made  in 
Section  3.8  that  the  dominant  term  in  the  expression 
for  round-off  error  using  Schur's  algorithms  is  of 
higher  order  than  the  corresponding  term  for  Gaussian 
elimination . 

3.  The  predicted  and  actual  values  of  the  round¬ 
off  errors  incurred  in  computing  the  inverse  of  a  block- 
symmetric  matrix  and  solution  of  a  system  of  equations 
using  Charmonman ' s  algorithms  are  smaller  compared  to 
those  obtained  by  using  Schur's  algorithms  and  Gaussian 
elimination ,  in  most  cases.  This  suggests  that  Charmonman’ s 
algorithms  are  less  susceptible  to  the  accumulation  of 
round-off.  These  remarks  are  in  complete  agreement  with 
the  observations  made  in  Sections  4.5  and  4.8  that  the 
dominant  term  in  the  expression  for  round-off  error 

using  Schur's  algorithms  is  of  higher  order  than  the 
corresponding  term  for  Charmonman ’ s  algorithms. 

4.  The  predicted  and  actual  round-off  errors  in 
inverting  band  matrices  and  solving  systems  of  equations 
with  band  matrices  are  almost  equal,  whereas  the  predicted 
round-off  error  using  Gaussian  elimination  without  minimum 
storage  is  comparatively  larger  than  the  actual  value. 

This  is  due  to  the  fact  tnat  only  a  few  arithmetic  operations 
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are  Involved  in  Inverting  band  matrices  and  solving 
systems  of  equations  with  band  matrices,  and  that  all 
computations  using  APL\360  are  done  in  double-preci¬ 
sion  . 

7 . 2  Suggestions  for  Future  Research 

1.  The  error  bounds  obtained  in  Chapters  III,  IV 
and  V  could  be  significantly  improved  if  the  analysis  is 
carried  out  using  accumulation  of  inner  products  in 
floating-point, 

2.  Winograd  [20],  has  proposed  a  new  algorithm 
for  matrix  multiplication  which  requires  smaller  number 
of  multiplications .  This  algorithm  can  be  applied  to 
invert  matrices  and  solve  systems  of  equations  using 
methods  described  In  Chapters  III,  IV  and  V,  provided 
partitioned  matrices  and/or  whole  matrices  are  inverted 
using  a  variation  of  Gaussian  elimination  given  In 

[3,  pp .  136-137]  or  [20].  Perhaps,  better  error  bounds 
for  the  predicted  round-off  error  Incurred  in  the 
computation  of  inverse  of  a  matrix  and  the  solution  of  a 
system  of  equations,  could  be  obtained  using  the  new 
algorithm  for  inner  products. 
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APPENDIX 

LISTINGS  OP  FORTRAN  IV 
AND  APL\360  PROGRAMS 

Listings  of  FORTRAN  IV  subprograms  and  APL\360 
functions  for  inverting  matrices  and  solving  systems  of 
linear  algebraic  equations  are  included  in  the  following 
pages.  The  SUBROUTINE  GAUSSC  computes  the  inverse  of  a 
matrix  or  solution  of  a  system  of  linear  equations  using 
Gaussian  elimination  with  complete  pivoting.  The  SUBROUTINE 
ERRBD2  j  SUBROUTINE  ERRBD4  and  SUBROUTINE  ERRBD5  compute 
the  inverse  of  a  general  and  block-symmetric  matrix.  The 
SUBROUTINE  ERRBD3  and  ERRBD6  compute  the  solution  of  system 
of  equations  with  general  and  block-symmetric  coefficient 
matrices.  The  SUBROUTINE  SCHURI  constructs  the  inverse  of 
a  partitioned  matrix.  The  SUBROUTINE  MADD  and  SUBROUTINE 
MPROD  compute  the  sum  and  product  of  two  matrices >  respectively. 
The  SUBROUTINE  NORM  computes  the  °°-norm  of  an  arbitrary 
matrix  and/or  the  element  of  maximum  modulus  of  a  vector. 

The  APL  function  BAND  computes  the  inverse  and  solution 
of  equations  using  minimum  storage  for  band  matrices. 

The  unit  round-off  error  for  IBM  System  360/67 

-21  -53 

computer  is  2  for  single-precision  arithmetic  and  2 
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for  double-precision  arithmetic.  Although  the  word  length 
in  single  and  double-precision  arithmetic  is  32  and  6H 
bits,  1  bit  is  used  for  the  sign,  7  bits  for  the  exponent, 
and  three  leading  bits  of  the  fractional  part  may  be  zero. 


since  the  normalization  is  done  in  hexadecimal. 
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SUBROUTINE  GAUSSC l N,M , ArX , SUL,G) 

DIMENSION  A l 20,20)  ,X( 20,20} , SQL (20,20)  , NCI  20) , NR ( 20) 
C 

C  THIS  SUBROUTINE  COMPUTES  THE  INVERSE  OE  A  MATRIX  OR 
C  SOLUTION  OF  A  SYSTEM  OF  EQUATIONS  CF  ORDER  N  BY 
C  GAUSSIAN  ELIMINATION  WITH  COMPLETE  PIVOTING. THE  INVERSE 
C  OR  SOLUTION  VECTOR  IS  RETURNED  VIA  THE  ARRAY  SOL. 

C  M  IS  THE  NUMBER  CF  RIGHT-HAND  SIDES. 

C  X  IS  THE  MATRIX  OF  RIGHT-HAND  SIDES. 

CGIS  THE  ELEMENT  OF  MAXIMUM  MAGNITUDE  AT  ALL  STAGES  IN 
C  THE  REDUCTION  OF  A  INTO  UPPER-TRIANGULAR  MATRIX. IT  IS 
C  USED  IN  THE  COMPUTATION  OF  THE  PREDICTED  RELATIVE  ERROR 
C  FOR  VARIOUS  ALGORITHMS. 

C 

COMMON  FLAG 
REAL  MU 

I F { N . EQ • I ) GO  TO  35 
DECLARE  ZERO  CRITERION 
ZERQ=lE-8 

CHECK  THE  NUMBER  OF  EQUATIONS  GIVEN. 

IFIM.EQ.DGO  TO  4 
DO  2  I = 1 , N 
DO  2  J=1,N 

2  XfI,J)=0 
DO  3  I  =  1 » N 

3  xu,n=i 

4  NM 1=N- I 
NP1=N+1 

THE  BEGINNING  OF  FORWARD  ELIMINATION. 

DO  31  1=1, NM1 
IP  1=1  +  1 
NR(i)=I 
NC ( I )=I 

USE  OF  COMPLETE  PIVOTING. 

PIVOT=C 
DO  5  J=I,N 
Du  5  K=  I , N 

I F ( PIVOT. GE.ABSI A( J,K)  )  JGu  TG  5 
P I VOT=  ABS ( A { J  »  K )  ) 


THE  VECTORS  NR  AND  NC  NOTE  THE 
I  NT  E  RCH ANG  E$,R  ESP  EC  T I  V ELY ♦ 


NECESSARY  ROW  AND  COLUMN 


' 


' 


v 
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NRt I ) = J 
NC (  I  I  =  K 
5  CONTINUE 

IF, AT  ANY  STAGE  OF  REDUCTION  OF  A  ,A  PIVOTAL  ELEMENT  IS 
LESS  THAN  THE  ZERO  CRITERION  , THE  MATRIX  IS  FOUND  TO  BE 
SINGULAR. 

IFt P  IVOT.LE. ZEROIGO  TO  65 

INTERCHANGE  ROWS, IF  NECESSARY. 

I F { N R ( I ) . EG. I ) GO  TO  15 
DO  10  J=1,N 
T  E  MP  =  A { I , J ) 

At  I  ,  J)=AtNRt I  )  ,  J  I 
A ( NR ( I )  ,  J)=TEMP 

10  CONTINUE 

DO  11  J=1,M 
T  EMP  =  X ( I ,J) 

X(I,J)=X{NR(I),J) 

X ( NR  t I ) , J)=TEMP 

11  CONTINUE 

INTERCHANGE  COLUMNS  ,IF  NECESSARY. 

15  I  F { N C  t  I  )  .  E Q . I ) GO  TO  25 
DO  20  J=1,N 
TEMP=A( JfI ) 

At J, I)=A(JfNC( I)) 

At  J, NC( I ) )=TEMP 
20  CONTINUE 
25  DO  31  J=IP1,N 

MIJ=-A{ J, I ) / A  (  I, II 
DO  30  K=1P1,N 
A! J,K}=A( J,K)+MIJ*A(I ,K) 

30  CONTINUE 

DO  31  L=1 »  M 

X{  J,U  =  X(  J,U +MI  J*X(I  ,  L  I 
1  CONTINUE 

IF, THE  A(N,N)  ELEMENT  OF  A  AFTER  N-l  REDUCTIONS  IS  LESS 
THAN  ZERO  CR I T ER I GN, THE  MATRIX  A  IS  FOUND  TO  BE 
SINGULAR. 

5  IF? A8S( A(NtN) I .LE.ZEROIGO  TO  70 

THE  BEGINNING  OF  BACK-SUBSTITUTION. 

DO  40  1=1, M 
X  t  N  »  I )  =  X  {  N ,  I  )  /  A  t  N » N  ) 


■ 


. 


' 
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40  CONTINUE 

I F ( N . E Q . 1 ) GO  TO  75 
DO  5  0  l<=2,  N 
I =NP 1-K 
I  P  1=  I  1 
DO  50  L=  1 ,  M 
SUM=0 

DO  45  J= I P 1 »  N 
SUM=SUM+A (  I,  J)*X(J*U 
45  CONTINUE 

X( I ,L)  =  (X(  I  »  L  )  -SUM ) / A ( 1,1) 

50  CONTINUE 

BACK- INTERCHANGE  OF  COLUMNS. 

DO  60  K= 1 , NM 1 
I  =N-K 

IF { NCI  I  )  •  EQ . I ) GO  TO  60  , 

DO  55  J=  1 »  M 
TEMP  =  X ( I , J) 

xi i , j)=x(ncii ) , j  ) 

X( NC( I ) , J)=TEMP 
55  CONTINUE 

60  CONTINUE 

STORE  INVERSE  OR  SOLUTION  VECTOR  IN  ARRAY  SOL. 

DO  61  1=1, N 
DO  61  J= 1 , M 

61  SOL  I  I , J )  =  X ( 3  » J } 

COMPUTE  ELEMENT  OF  MAXIMUM  MODULUS  IN  THE  REDUCTION  OF  A. 
G=0 

DO  62  1=1, N 
DO  6  2  J=  I  »  N 

62  G= AM AX 1(G,ABS(A(I,J))) 

GO  TO  75 

65  FL AG=1 

GO  TO  75 
7C  FL AG=2 

75  RETURN 
END 


■  ■ 


‘ 


■ 
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SUBROUTINE  ERRB01 { N,NRA» A, B,C, D,  RE NORM  ,  TNORM , HC  » 

L  PROD , EN2 ) 

DIMENSION  A (20,20) ,B{  20,20 ),C( 20,20) ,D( 20,20), 

1HC(  2  0,20)  ,X(  20,20  ,  AC  I  (20,20)  ,  PROD  (20, 20)  , 

2 PROD  1  ( 20,20) , DELTA (20 ,20) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  MATRIX  HC  USING  SUBMATRICES 
C  A , 8 , C  AND  D  AND  COMPUTES  AN  UPPER  BOUND  FOR  THE 
C  ABSOLUTE  ROUND-OFF  ERROR  EN2  INCURRED  IN  THE 
C  CALCULATION  OF  HC.THE  COMPUTED  HC  AND  EN2  ARE  RETURNED 
C  VIA  HC  AND  EN2.THE  PRODUCT  AC  I  $  6 ,  ‘WHERE  ACI  IS  THE 
C  COMPUTED  INVERSE  OF  A  IS  ALSO  RETURNED  VIA  ARRAY  PROD 
C  FOR  FUTURE  USE. 

C  N  IS  THE  ORDER  OF  THE  MATRIX  R. 

C  NR A  IS  THE  ORDER  Or  THE  MATRIX  A. 

C  RENORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  R. 

C  TNORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  A. 

C  EP  IS  THE  UNIT  ROUND-OFF  ER-ROR  FOR  IBM  360/67  COMPUTER. 

C 

DOUBLE  PRECISION  EP 
EP  =  2  .D0**(-21 ) 

NMR-N-NRA 


BNORM , CNCRM  AND  DNQRM  IS  THE  INFINITY  NORM  OF  MATRICES 
B,C  AND  D, RESPECTIVELY . 

CALL  NORM! NRA,NMR,8, BNORM) 

CALL  NCRM( NMR , NR A , C , CNGRM ) 

CALL  NORM( NMR , NMR, D, DNORM) 

COMPUTE  ACI, THE  INVERSE  OF  A. 

CALL  GAUSSC ( NRA , NR A , A, X , AC  I ,  G) 

ACNORM  IS  THE  INFINITY  NORM  OF  ACI. 

CALL  NGRM( NRA, NRA, ACI , ACNORM) 

EN1 1^2 .005*NRA**2+NRA**3 
EN12=EN11*TN0RM*G 

EN13=EN12+1 .+2.*NRA+NRA*INRA+2. ) *EP 

EN14=EN13*CNGRM*ACNQRM 

EN1=E?*(DN0RM+EN14*BN0RM) 

COMPUTE  PROD ,HC  AND  EN2 • 

CALL  MPROOI NRA, NRA, NMR, ACI , 8, PROD) 

CALL  MPROD ( NMR , NRA , NMR , C , PROD , PROD  1 ) 

CALL  M ADD l NMR , NMR , D , PR0D1 ,DELT  A) 

CALL  GAUSSC<NMR,NMR,DELTA,X,HC,G) 

CALL  NORM( NMR ,NMR , HC , HCNORM ) 
EN21=2.005*NMR**2+NMR**3 


' 


' 

■ 

‘ 


■ 


. 


138 


EN22=EP*EN21*G*HCNORM 
EN23  =  RENGRM:!!EN1 

EN2=(EN22+EN23)*lRENURM/( 1 . - RENORM*EN 1 ) ) 

RETURN 

END 


. 


1 

■ 


' 


. 
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SUBROUTINE  ERR8D2 l N , NR A ,  A  ,  B,  C  ,  D, RE NORM , TNORM , RC I , 

1  ERRS ) 

DIMENSION  A (20,20) ,8 ( 20,20) ,C( 20,20) ,0(20,20) , 

IRC  I (  20,20) , PROD (20, 20)  ,PR0D1 ( 20 , 20 )  , EC ( 20 , 20  )  , 
2FC(20,20) ,GC( 20,20)  ,HC  (20,20) , AC  I ( 20 ,20) ,X( 20,20) , 

3  S ( 20,20) ,DC( 20,20) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  INVERSE  OF  R  USING  SCHUR'S 
C  ALGORITHM  AND  CALCULATES  UPPER  BOUND  FOR  THE  ROUNU-OFF 
C  ERROR  INCURRED  IN  THE  C ALCULA T I  ON . THE  COMPUTED  INVERSE 
C  IS  RETURNED  VIA  THE  ARRAY  RC I • THE  PREDICTED  RELATIVE 
C  ERROR  IS  STORED  IN  ERRS. 

C  N  IS  THE  ORDER  OF  THE  MATRIX  R. 

C  NR A  IS  THE  ORDER  OF  MATRIX  A. 

C  A , B , C  AND  D  ARE  SU3MATRICES  OF  R. 

C  RENORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  R. 

C  TNORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  A. 

C  EP  IS  THE  UNIT  ROUND-OFF  ERROR  FOR  IBM  360/67  COMPUTER. 

C 

DOUBLE  PRECISION  EP 
EP=2.DC**(-21 ) 

NMR=N-NRA 
NR AP  1  =  NR A+l 
DO  5  1=1, NR A 
DO  5  J  =  1 »  NRA 
5  S( I , J)  =  A( I  ,J) 

STORE  A  FOR  FUTURE  USE. 


COMPUTE  HC. 


CALL  ERRBD1 (N, NRA, A, B,C,D, RENORM, TNORM , HC , PROD, EN2 ) 


G  IS  THE  ELEMENT  OF  MAXIMUM  MODULUS  IN  THE  REDUCTION 
OF  A. 


CALL  GAUSSC( NRA, NRA, S, X, AC  I  ,  G) 
CALL  NORM (NRA, NR A, AC  I , AC NORM) 
CALL  NQRM( NRA,NMR,6»3NGRM) 

CALL  NCRMC NMR , NRA , C , CNGRM ) 


HCNGRM  IS  THE  INFINITY  NORM  OF  THE  MATRIX  HC . 


CALL  NGRM{ NMR , NMR , HC , HCNQRM ) 

EN11=2.005*NRA**2+NRA**3 

EN12=EN11*TN0RM*G 

EN31=TN0RM*BNGRM*EN2 

EN  32= ( EN12+N  +  NRA-NMR) *EP 

EN3=EN31+EN32*ACNORM*BNQRM*HCNORM 

EN41=EN2+EN12*RENQRM+EP* (N*NRA*NMR*EP ) *HCNORM 


' 


' 


" 

' 

' 
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' 
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EN4=EN41*ACN0RM*CN0RM 


COMPUTE  GC. 

CALL  MPRGD(NM.R,NMR,NRA,HC,C, PR0D1) 

CALL  MPROD ( NMR , NRA » NRA »  PR0D1 ,ACI ,GC) 

DO  10  1=1 » NMR 
00  10  J  =  l ,  N  R  A 
10  GC( I » J )=-l . *GC (  I  i J) 

GC  NORM  IS  THE  INFINITY  NORM  OF  THE  MATRIX  GC 

CALL  N0RM( NMR,NRA, GC, GCNCRM) 

EN51=TN0RM*BNGRM*EN4 

E  N  5  2  =  {  1 .  +  EM  2 )  *ACNORM 

EN53=N*(  1  • *EP )  +NRA*NMR*EP 

EN54= ( EN52+EN53*ACNQRM) *BN0RM*GCN0RM 

EN5=EP*(EN5L+EN52  +  EN54). 

COMPUTE  RELATIVE  ERROR  ERRS. 


ERR  1  =  EN2  +  EN4 
ERR2=EN3+EN5 
ERRS  =  A  MAXI ( ERR1, ERR2) 

ERRS= ERRS/ RENO RM 

COMPUTE  FC 

CALL  MPROD (NRA, NMR, NMR, PROD, HC ,FC) 

DO  20  1=1, NRA 
DO  20  J=1 , NMR 
20  FC( I , J)=-1.*FC( I , J) 

CGMPUTE  EC. 

CALL  M PROD ( NRA, NMR, NRA, PROD, GC , PR001 ) 
CALL  MADDl NRA, NRA, AC  I , PR0D1,  EC) 

DO  103  1=1, NMR 
DO  103  J  =  1 » NMR 
DC (  I  ,J)  =  HC(I,J) 

103  CONTINUE 

COMPUTE  RCI . 

CALL  SCHURK  N  ,  NRA  ,  EC ,  F  C  ,GC ,  DC ,  RC  1 1 

RETURN 

END 


■ 

. 
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SUBROUTINE  ERRB03 ( N , NR A , A , B , C »  0 , RHS , TN QRM , RENORM, 
1XEMAX,XC1 ,ERRS) 

DIMENSION  A  (20,20  ,  3(  20, 20  )  ,  C  (  20 , 20  )  ,  D  (  20 , 20  )  , 

1RHS( 20,20 ,XC2(20,20) ,X(20,20) ,ACI (20, 20) ,B1 (20,20) , 
292(23, 20), PROD  120, 20), PROD 1(20,20)  , PR0D2 ( 20 , 20 ) , 
3DELTA( 20,20) , DEL (20,20) ,XC1( 20,20) ,001(20,20), 
4RS(2Q,20),S(20,20) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  SOLUTION  OF  A  SYSTEM  OF 
C  EQUATIONS  USING  SCHUR'S  ALGORITHM  AND  CALCULATES  AN 
C  UPPER  BOUND  FOR  THE  ROUND-OFF  ERROR  INCURRED  IN  THE 
C  CALCULATION  .THE  COMPUTED  SOLUTION  IS  RETURNED  VIA  ARRAY 
C  XC1  AND  THE  PREDICTED  RELATIVE  ERROR  IS  STORED  IN  ERRS. 

C  N  IS  THE  ORDER  OF  THE  COEFFICIENT  MATRIX  R. 

C  NR A  IS  THE  ORDER  OF  THE  MATRIX  A. 

C  A,B,C  AND  D  ARE  THE  SUBMATRICES  OF  R. 

C  RHS  IS  THE  MATRIX  OF  RIGHT-HAND  SIDES. 

C  TNORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  A. 

C  RENORM  IS  THE  INFINITY  NORM  UF  THE  EXACT  INVERSE  OF  R. 

C  XEMAX  IS  THE  INFINITY  NORM  OF  THE  EXACT  SOLUTION  OF  THE 
C  GIVEN  SYSTEM  OF  EQUATIONS. 

C  EP  IS  THE  UNIT  ROUND-OFF  ERROR  FOR  THE  IBM  360/67 
C  COMPUTER. 

C 

DOUBLE  PRECISION  EP 
EP=2.D0**(  -21  ) 

NMR=N-NRA 


STORE  THE  MATRIX  A  FOR  FUTURE  USE. 

DO  5  1=1, NRA 
DO  5  J= 1 , NRA 
S( I , J)=A( I  , J) 

5  CONTINUE 

ANORM  IS  THE  INFINITY  NORM  OF  THE  MATRIX  A. 

CALL  NORM! NRA, NRA, A, ANORM) 

PARTITION  THE  RIGHT-HAND  SIDE  RHS. 

DO  10  1=1, NRA 
10  B 1 ( 1 , 1 )  =  RH  S ( 1 , 1 ) 

NRAP 1=NRA+ 1 
DO  20  I=NRAP1,N 
1 1  =  1 -NRA 

20  B2( I1,1)  =  RHS(  1, 1) 

AC  I  IS  THE  COMPUTED  INVERSE  OF  A. 

CALL  GAUSSC ( NRA, NRA, A»X, AC I ,G) 


' 

‘ 
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CALL  NORM { NRA , 1 , B1  ,  AMAX61 ) 

CALL  NORM! NMR , 1 , B2  , AM AX 62 ) 

EN 11=2 .0C5*NRA**2+NRA**3 
EN12-EN11*TNGRM*G 

EN13=EN12+1 *  +2 • *NRA+NR A* ( NR A+2 . ) *EP 

CNGRiM  IS  THE  INFINITY  NORM  OF  THE  MATRIX  C. 

CALL  NGRMl NMR , NRA , C , CNCRM) 

CALL  NORM (  NR  A , NR  A , ACI , ACNORM) 

EN  14=EN1  3*  C  NORM5*  ACNORM 
ENS2=EP*I AMAXB2+ EN 14* AMAXB 1 ) 

COMPUTE  THE  COEFFICIENT  MATRIX  AND  THE  RIGHT-HAND  SIDE 
FOR  COMPUTING  THE  SOLUTION  VECTOR  OF  ORDER  N-NRA. 

CALL  MPROD ( NMR , NRA , NR A , C , AC  I , P ROD ) 

CALL  MPROD (NMR, NR A, 1 , PROD, 61 , PROD1 ) 

CALL  MADD( NMR, 1, 62, PRO 01, X) 

CALL  MPROD (NMR, NRA, NMR, PROD, 6, PR0D1) 

CALL  M ADD ( NMR NMR, D,PRGD1, DELTA) 

CALL  GAUSSC ( NMR , 1 , CELT A , X , XC 2 , G ) 

CALL  NORM( NMR, 1 »  XC2 , AM AXC2 ) 

EK1=EP*EN1 1*G 

DNORM  IS  THE  INFINITY  NORM  OF  MATRIX  D. 

CALL  NORM l NM R , NM R , D , ON ORM ) 

EN 1=  E  P* ( ON OR M+  EN 1 4*BN0 RM ) 

EN3= (ENS2-M  EN1+EK1 ) * AM AXC2 ) *RENQRM 
I  F ( NR A-NMR ) 15,21, 15 

COMPUTE  THE  SOLUTION  VECTOR  OF  ORDER  NRA, IF  NRA  NOT 
EQUAL  TO  N/2. 

15  EN41  =  1 .  +  2  .  *:NMR  +  NMR*  (NMR  +  2  # )  *EP 

DENORM=RENORM 

CALL  GAUSS C( NMR, NMR, D, X,DCI ,G) 

CALL  NORM (NMR, NMR, DC  I , DCNORM ) 

EN42-2.C05*NMR**2+NMR**3 

EN43  =  (  G*EN42*DEN0Rrt+EN41 }  *BNQRM*DCNQRM 

EN4= EP* ( ANQRM+EN43-CN0RM) 

EN5=EP*( AMAXB1+EN43*AMAXB2) 

CALL  MPROD (NRA, NMR, NMR , B , DC  1 ,PRGD) 

CALL  MPROD (NR  A, NMR, NRA, PROD, C, PROD  I) 

CALL  MPROD (NR A, NMR , 1 , P ROD , B2 , PR0D2 ) 

CALL  M A 0 D ( NRA ,NR A , S , PR0D1 , DEL ) 

CALL  MAOD ( NRA, 1, Bl, PRCD2,  RS  > 

COMPUTE  THE  SOLUTION  VECTOR  OF  ORDER  NRA. 


' 
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CALL  GAUSS C ( NR  A  1 1 , DEL , RS , XC 1 , G ) 

CALL  NORMCNRA, 1 , X Cl , AM AXC1 ) 

EK2=EP*EN1 1*G 

E  N6= ( EN5+ ( EN4+ER2 ) *AMAXC1 ) *RENGRM 
GG  TO  25 

COMPUTE  THE  SOLUTION  VECTOR  OF  ORDER  N/2.IF  NR A  EQUAL 
TO  N/2. 

21  CENORM=RENORM 

CALL  MPRODCNMRtNMR, 1,D»XC2, PROD) 

CALL  MADD(NMR,1» B2»PROD,RS) 

CALL  GAUSS C ( NMR  *  1 »C»RS»XC1»G) 

CALL  NORM! NMR, 1 , XC 1 , AM AXC 1 ) 

EN6= ( EN4+EK2*AMAXC1 ) *C  ENORM 
25  ERRS=AMAX1 ( EN3  »  EN6  ) 

COMPUTE  THE  PREDICTED  RELATIVE  ERROR  ERRS. 

ERRS=ERRS/XEMAX 

CONSTRUCT  SOLUTION  VECTOR  OF  ORDER  N. 

DO  30  I  =  NR  A  P 1 » N 
XC1( I ,  1)=XC2(  I-NRA, 1 ) 

RETURN 
END 


30 
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SUBROUTINE  ERRBD4 ( N , A , B ,C , D , RE NORM , RC I ,  ERRS ) 
DIMENSION  A (20,20) ,  B( 20,20 ) , RC 1(20,20) ,EC( 20,20) , 

1 FC ( 20,20) , PROD (20, 20) , S(20, 20 ),C(20, 20 ),D{ 20, 20), 

2HC (20,20) 

C 

C  THIS  SUBROUTINE  COMPUTES  THE  INVERSE  OF  4  BLUCK- 
C  SYMMETRIC  MATRIX  USING  SCHUR'S  ALGORITHM  AND  CALCULATES 
C  AN  UPPER  BOUND  FOR  THE  ROUND-OFF  ERROR  INCURRED  IN  THE 
C  CALCULAT  ION. THE  INVERSE  OF  THE  BLOCK-SYMMETRIC  MATRIX 
C  IS  RETURNED  VIA  THE  ARRAY  RCI  AND  THE  PREDICTED  RELATIVE 
C  ERROR  IS  STORED  IN  ERRS. 

C  N  IS  THE  ORDER  OF  THE  MATRIX  R. 

C  A , B , C  AND  D  ARE  SUBMATRICES  OF  R. 

C  RENORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  R. 

C  EP  IS  THE  UNIT  ROUND-OFF  ERROR  FOR  THE  IBM  360/67 
C  COMPUTER. 

C 

DOUBLE  PRECISION  EP 
EP  =  2  .D0^*( -21 ) 

NR A=N/2 
NMR=  N/2 
NRAP  L=NR A+ 1 

STORE  MATRIX  A  FOR  FUTURE  USE. 

DO  5  1=1, NRA 
DO  3  J= 1 » NRA 
S(  I  ,  J  )  =  A C  I  ,  J  ) 

5  CONTINUE 

TNORM=RENORM 


COMPUTE  MATRIX  EC. 

CALL  ERR6D1 (N ,NR A, A , 3 , C , D, RENORM , TNORM , HC , PROD ,  EN2 ) 
DO  10  1=1, NMR 
DO  10  J=1 , NMR 
10  EC (  I  ,  J  )  =HC ( I , J) 

E2=EN2 

COMPUTE  MATRIX  FC . 

CALL  ERR9D1 ( N ,NR A, D, C , B , S, RE NORM, TNORM ,HC , PROD,  EN2 ) 
E4=EN2 

DO  20  1=1, NMR 
DO  20  J=1 , NMR 
20  FC ( I  ,  J  )  =HC ( I , J ) 

ERRS=E2+E4 

COMPUTE  RELATIVE  ERROR  ERRS. 


ERRS=ERRS/ RENORM 


' 
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CONSTRUCT  INVERSE  RC I  USING  EC  AND  FC. 

CALL  SC MUR I ( N,NRA, EC, FC, FC, EC , RC I ) 

RETURN 

END 


- 


. 
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SU8RCUT  I  FIE  ERR6D5(N,A,a,REN0RM,RC,ERRC) 

DIMENSION  A (20,20) ,U( 20,20 ),RC( 20,20) ,P (20, 20), 
1Q(20,20),PCI (20,20)  ,QCI  (20,20  ,X(20,20) , EC (2 0,20), 
2FC ( 20,20) 

C 

C  THIS  SUBROUT  E.NE  COMPUTES  THE  INVERSE  OF  A  BLOCK- 
C  SYMMETRIC  MATRIX  USING  CHA  RMUNMAN*  S  ALGORITHM  AND 
C  CALCULATES  A;l  UPPER  BOUND  FUR  THE  ROUND-OFF  ERROR 
C  INCURRED  IN  THE  C A LC ULA T I  ON . THE  INVERSE  IS  RETURNED  VIA 
C  THE  ARRAY  RC  AND  THE  PREDICTED  RELATIVE  ERROR  IS  STORED 
C  IN  ERRC. 

C  N  IS  THE  ORDER  OF  THE  BLOCK-SYMMETRIC  MATRIX. 

C  A  AND  B  ARE  THE  SUBMATRICES  OF  R  OF  ORDER  N/2  EACH. 

C  RENORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  R. 

C  EP  IS  THE  UNIT  ROUND-OFF  ERROR  FOR  THE  IBM  360/67 
C  COMPUTER. 

C 

DOUBLE  PRECISION  EP 
EP=2.D0**C -21 ) 

NRA=N/2 

CALL  NQRHl NR  A , NR A, A.ANORM) 

CALL  NORM NRA.NRA, B,BNORM) 

EN1=EP*( ANCRM+3NGRM) 

EN2=  2 . *E P* ( ANORM+2 . *BNQRM ) 

EN31=2.005*NRA**2+NRA**3 

COMPUTE  P  AND  Q. 

DO  10  1=1, NRA 
DO  10  J=  1 , NRA 
P(I , J ) = A { I , J ) +B (  I, J) 

10  Q(I»J)=P(I»J}-2.*B(I»J) 

PCI  IS  THE  COMPUTED  INVERSE  OF  P. 

CALL  GAUSSCt  NRA,NRA,P, X,PCI ,  G) 

CALL  NORM (NRA, NRAtPCI, PCNORM) 

PEN0RM=2.*REN0RM 
GENORM=2.*RENORM 
EN32=G*EN3 1*PCN0RM+PEN0RM*EN 1 
EN3- EP*EN323M PENORM/ ( 1 . -PENO RM*EN 1 ) ) 

QCI  IS  THE  COMPUTED  INVERSE  OF  Q. 

CALL  GAUSS C ( NR A, NR A, Q, X, QC I , G) 

CALL  NORM! NRA, NRA, QCI , QCNORM ) 
EN41=G*EN31*GCN0RM+CENGRM*EN2 
EN4=EP*EN41*( QENORM/ ( 1 • -UEN0RM*EN2 ) ) 

EN  3=0 . 5* { EN3+EN4 HO. 5*EP*t PCNORM+QCNURM) 

ERRC=EN4+2 .*EN5 
C 


' 


' 


. 


,  } 
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COMPUTE  THE  PREDICTED  RELATIVE  ERRCR  ERRC. 

ERRC=ERRC/RENCRM 

COMPUTE  EC. 

DO  20  1=1, NRA 
DO  20  J=1,NRA 

20  EC  (  I  ,  JJ  =0,5*  (PC  I  (  I  ,  J)+GCI  (  I,  J)  } 

COMPUTE  FC. 

CALL  MADOI NR A , MR A , EC , QC I ,FC) 

CONSTRUCT  INVERSE  OF  R  USING  EC  AND  FC . 

CALL  SCHURKN,  NR  A ,  EC ,  FC  ,  FC ,  EC,RC) 

RETURN 

END 


) 


■ 


' 


■ 


. 
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SUBROUTINE  ERR6D6  {  N  ,  A  >  B,RHS,  RENORM ,XEMAX , XC2 , ERRC ) 
DIME  NS IOW  A ( 20, 20 )  ,6(20,20) , RMS (20, 20)  ,YC2( 20,20), 
1B1( 20,201, B2 l 20,20) , RSI (20,20 ) ,RS2( 20,20 ) ,P( 20,20) , 
2U<20,20),YC1( 20,20) ,XC 1(20,20)  ,XC2( 20, 20)  ,SCL( 20,20) 
C 

C  THIS  SUBROUTINE  COMPUTES  THE  SOLUTION  OF  A  SYSTEM  OF 
C  EQUATIONS  WITH  BLOCK-SYMMETRIC  COEFFICIENT  MATRIX  USING 
C  CHARMONMAN *  S  ALGORITHM  AND  CALCULATES  AN  UPPER  BOUND 
C  FOR  THE  ROUND-OFF  ERROR  INCURRED  IN  THE  CALCuLAT ION . THE 
C  COMPUTED  SOLUTION  IS  RETURNED  VIA  THE  ARRAY  XC2  AND  THE 
C  PREDICTED  RELATIVE  ERROR  IS  STGkED  IN  ERRC • 

C  N  IS  THE  ORDER  OF  THE  COEFFICIENT  MATRIX  R. 

C  A  AND  6  ARE  THE  SUBMATRICES  OF  R. 

C  RHS  IS  THE  MATRIX  OF  RIGHT-HAND  SIDES. 

C  RENORM  IS  THE  INFINITY  NORM  OF  THE  EXACT  INVERSE  OF  R. 

C  XEMAX  IS  THE  INFINITY  NORM  OF  THE  EXACT  SOLUTION  OF  THE 
C  GIVEN  SYSTEM  OF  EQUATIONS-. - 

C  EP  IS  THE  UNIT  ROUND-OFF  ERROR  FOR  THE  IBM  360/67 
C  COMPUTER. 

C 

DOUBLE  PRECISION  EP 
EP=2.D0**{-21) 

NRA=  N/2 
NR AP 1=NRA+ 1 


PARTITION  RIGHT-HAND  SIDE  RHS. 

DO  10  1  =  1 , NRA 
B1U,1  )=RHS  ( l  ,1) 

I 1=NRA+I 

10  B2(  I  ,li=RHS( 11,1 ) 

DO  20  1=1, NRA 

RSI (  1 , 1 )  =  B1(  I  ,  1) +B2( I , 1) 

20  RS2(  I,1)=RS1( 1,1 )-2. *82(1,1) 

COMPUTE  P  AND  Q. 

DO  30  1=1, NRA 
DO  30  J  =  1 »  NRA 
P ( I , J)=A(I ,J)+3( I,J) 

30  Q( I , J)=P( I ,J)-2.*B(I , J  ) 

CALL  NORIK  NRA, NRA, A, ANORM) 

CALL  NORM  NR A , NR A, S , BNORM ) 

PEN0RM=2.*RENCRM 
QENQRM=2.*RENQRM 
CALL  NORM NRA, 1,B1,AMAX31) 

CALL  NORM ( NRA, 1 , 52  »  AMAXB2 ) 

COMPUTE  SOLUTION  OF  THE  SYSTEM  OF  EQUATIONS  WITH 
COEFFICIENT  MATRIX  P. 


, 

■ 

* 
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' 


■ 


: 


' 
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CALL  GAUSSC(NRA,1,P,RS1,YC1,G> 

CALL  NQRM( NRA, 1, YC1 iAMAXYI) 

EN51=(  G*{  2.0Q5*NRA**2+NRA**3)  + ANGRM+BNGRM ) *AMAXY  1 

EN52=AMAX81+AMAXB2 

EN5=EP*(EN51+EN52)*PEN0RM 

COMPUTE  SOLUTION  OF  THE  SYSTEM  OF  EQUATIONS  WITH 
COEFFICIENT  MATRIX  Q. 

CALL  GAUSSCINRA, 1,Q,RS2,YC2, G) 

COMPUTE  FIRST  N/2  COMPONENTS  OF  THE  SOLUTION  VECTOR. 

DO  35  1=1 t NR A 

35  XC1( I,1)=0.5*(YC1{ I,1)+YC2<  Itl)) 

CALL  NORM!  NR  A  ,  1 ,  YC2  ,  A.M  AXY2  ) 

EN61=(G*( 2.005*NRA**2*NRA**3 )+2.*{ AN JRM+2 . *BNURMI ) 
1+AMAXY2 

EN6= EP*{ EN53  +  EN61  ) *QENGRM 

COMPUTE  LAST  N/2  COMPONENTS  OF  THE  SOLUTION  VECTOR. 

DO  40  1=1, NRA 

40  XC2(I,1)=XC1(  I,1)-YC2II,1) 

CALL  NORM ( NR A , 1 , XC 1 , AMAXD 1 ) 
ERRC=EN5+2.*EN6+EP*AMAXC1 

COMPUTE  THE  PREDICTED  RELATIVE  ERROR  ERRC. 

ERRC=ERRC/XEMAX 

CONSTRUCT  SOLUTION  VECTOR  OF  ORDER  N. 

DO  50  I  =  NR AP 1 , N 
50  XC2I  I,1)  =  XC1(  I-NRA, 1) 

RETURN 

END 


■ 


' 

v 
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SUBROUTINE  SCHURI l N , NRA t A , B , C , Q, RC I ) 

DIMENSION  A (20,20 ,B( 20,20 5 ,C(2C,20) ,0(20,20), 

1RCI  (  20,20) 

THIS  SUBROUTINE  CONSTRUCTS  A  MATRIX  OF  ORDER  N  FROM 

THE  SUBMATRICES  A,3,C  AND  D, AND  STORES  IT  IN  ARRAY  RCI. 

NRA  IS  THE  ORDER  OF  MATRIX  A. 

NRAP 1=NRA+ 1 
DO  10  1=1, NRA 
DO  10  J= 1 , NR A 
RCI  {  I , J )  =  A ( I , J) 

10  CONTINUE 

DO  15  1=1, NRA 
DO  15  J  =  NR AP 1 »  N 
RCI  (  I  »  J ) =3 ( I , J-NRA) 

15  CONTINUE 

00  20  I  =  NR  AP 1 , N 

DO  20  J= 1 , NR A 

RCI (  I  ,  J ) =C I  I -NRA, J) 

20  CONTINUE 

DO  25  I =NR AP 1 , N 
DC  25  J=NR AP 1 , N 
RCI  (  I  , J )  =  0( I-NRA, J-NRA) 

25  CONTINUE 
RETURN 
END 


■ 


' 


‘ 
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SUBROUTINE  NORM { N , M , A ,  XNGRM ) 

DIMENSION  A ( 20  »  20 ) 

THIS  SUBROUTINE  COMPUTES  THE  INFINITY  NORM  OF  AN 
ARBITRARY  MATRIX  A  OF  DIMENSION  N-BY-IM  AND  STORES  IT  IN 
XNCRM . 


XN0RM=0 
DO  20  1=1, N 
S  U  M=  0 

DO  10  J=1 , M 

10  SUM=  SUM  +  AB  S ( A ( I , J ) ) 

IF ( XNORM. GT • SUM) GO  TO  20 
XNORM=  SUM 
20  CONTINUE 
RETURN 
END 


, 


.  .  . 

‘ 


. 


- 


' 

' 
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SUBROUTINE  MPROD (N, M, K , A, B , AMLT ) 

DIMENSION  A (20, 20)  ,8( 20,20) , AMLT (20, 20) 

THIS  SUBROUTINE  COMPUTES  THE  PRODUCT  OF  TWO  MATRICES 
CONFORMABLE  FOR  MATRIX  MULT  I  PL IC AT  ION . THE  PRODUCT  IS 
RETURNED  VIA  THE  ARRAY  AMLT. 

DO  20  1=1, N 
DO  20  J=  I , K 
SUM=0 

DO  15  L  =  1 »  M 

15  SUM=SUM+A( I ,L)*B(L, J) 

20  AMLT ( I , J)=SUM 
RETURN 
END 


' 


i 


/ 
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SUBROUT  I NE  MADO ( N , M , A , B , ADD ) 

DIMENSION  A( 20,20) ,3(20,20) ,ADD( 20, 20) 

THIS  SUBROUTINE  SUBTRACTS  MATRIX  B  OF  DIMENSION  N-BY-M 
FROM  MATRIX  A  OF  THE  SAME  D I  MENS  I ON . TH E  RESULT  IS  STORED 
IN  ARRAY  ADD. 

DO  10  1=1, N 
DO  10  J- 1 ,  M 
ADD {  I , J)  =  A( I , J  } -8 {  I  ,  J) 

10  CONTINUE 
RETURN 
END 


s 

« 


■ 


' 


■ 
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V  BAND 

C  1  ]  -> (  (  0 <  (  ~ 2 x/,/ )  + 1  +  MC+ 1 +MUD  +MLD  )v(  ( 11 LD <  0  )  WUD <  0  )  )  / 

75 

[2]  MC+  L  111 , MC 

[  3  ]  MZ+(  ( MU+1  )  t  2  )  xMU+MC  -  ( MPP+ 1  ) 

[  4  ]  KA+UlxMC )  -  (  (/./£+ 1  )  v  2  )  xMLHIC- ( MLF+ 1  ) 

[5]  MR+M-ML 

[6]  ->(  ( MLD-  0  )  a  ( /■/£//)  =  0  )  )  /  6  2 

[7]  ->U1LD=  0)/47 

[8]  ->(/7Z  =  0)/13 

[93  I+-J+- 1 

[10]  [  i  MZ7Z7  +  J  ]  ,  (  ( MC -MUD + I )  p  0  )  ,  A  [  (  (  p  A  )  oj  (  p  A  )  -MUD+J  )  /  i  p  A] 

[11]  J+J+MC+ 1 

[12]  ->(MU>I+I+ 1)/10 

[13]  PIV<-[/\A 

[14]  TOL+EPSxPIV 

[15]  T+U1LD+  l)p  KL+KG+KI+Kl  0 

[16]  KM*- 1 

[17]  ■*(  I<M-MUD+1  ) /19 

[18]  KM+M-MUD 

[19]  +U4RZMLD+I)  /  2  5 

[20]  KL+KL+ 1 

[21]  +(KL<MUD)/ 25 

[22]  +  (  (pJP)<2)/47 

[23]  </«-!+  (  \  p  UK)  i  PP«-  (  u  [  JP  ]  )  x  PJ  7<-  r  /  u  [  JP«-J7-'  [  1+  x  (  p  c7P )  - 1  ]  ] 

[24]  ->-2  6 

[25]  «7<-J+(  xpJP)  xPP^-(  M[JP]  )  x  PIV+-  [  /  |/i[eTP+( -(T<-T+(  (MLP+l-KL) 
p  0  )  ,  X  KL  )  )  + 1  +/;px  (  ”  1  +  x  M )  [  KM+  X  PPZK  1  ]  ] 

[26]  ->( PIV<TOL  )  /  II 

[27]  -*(e7=I«-I+l  )/34 

[28]  TEMP*- A [ JJ+JK [ PP ] +  1 +  x  MC - P I ] 

[29]  .4  [  e7c7  ]  [  S'  -5-c7  P  [  1  ]  +  1  -i-  x  MC-KI  ] 

[30]  AIS^+TEMP 

[31]  TEMP+RLJil 

[32]  P  [  c7 ;  ]  *-R  [  J  ;  ] 

[33]  P[  J;  l*- TEMP 

[34]  K+ 2 

[35]  P[J;  >2?[I;  MAIJKUYS 

[36]  AlSl+AlS+JKl  1]+  l+xMC-KIItAUKllYl 

[37]  P[J+c7PxJP[P-l]  ;  >P[J+c7Pic/P[P-l]  il-AUKinixRLlil 

[38]  AlLLl+MAZLL+JKin+  1  +  x  M  C- -  K I  ]  -  A  [  JK  [  K  ]  ]  *  /,  [  J  K  [  1  ]  + 

1 + \MC-KI ] ) [ 1  +  x  MC-KI+ 1 ] ) , 0 

[39]  -*((pe7702J2)/43 

[40]  ->(  (pJK)<MLD)  /42 

[41]  ->(  (MLD+1  )>K*-K+1 )  /37 

[42]  ->(  (  (pJTO  +  l  )>R+K+1  )/37 


' 
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[43]  /  4  6 

[44]  KI+1CI+ 1 

[45]  ->(  (  (  p  JK  )  =//+ 1  -  J  )  A  (  J  <M-  1  )  )  /  2  2 

[46]  ->-(J</M)/16 

[47]  i?  [  ;  ]  +R  [ M ;  ]  iA  [  AM  -MLD  ] 

[48]  I  <-!!-! 

[49]  tJ-*-L*-X+-Q 

[50]  V*  i  1  +K*~  1 

[51]  i?[  J ;  >A?[ I ;  ]  -i4  [  II«-2  +  ( X- 1 )  +  +  /  (  ( N  9MCx  (  (  ( A/-  p <LD+  i  'T-’T, n 
+  1  )pl)  ,0  )  )[2+iA7-l]  )[J>  iM-1+Jl  ] x A?[ A:"-L  ;  ] 

[52]  F*-y,IT 

[53]  ^(L<AO/56 

[54]  K+K+l 

[55]  -+(X<L+L- 1)/51 

[56]  A? [  J  ;  ] «-i? [  J  ;  ]  t>1  [  7[  1  ]  -  1  ] 

[57] 

[58]  «W+1 

[59]  -*  (  ( M-  ( MLD  +  MUD  )  )  <  7>T  - 1 )  /  5  0 

[60]  ->(I<0)/65 

[61]  ->5  0,I>Z+1 

[62]  J«-l 

[63]  i?[  J ;  ]^i?  [  J ;  ]  [  J  ] 

[64]  +(M>I+I+ l)/63 

[65]  G-e-r/U 

[66]  EP+-2*  5  3 

[67]  F  l«-2  x  A7LP  x  (  3* MLD  )  +  (  3  x  AAt/C  )  +  1 

[68]  F2«-  ( AMD  + 1  )  x  U1LD+MUD  + 1 )  x  (  (  2x/;/LD  )  +/  047  +  4  ) 

[69]  F3x-0 . 5xFPx  (F1+F2  )x£ 

[70]  RCNORH+t /+/ \R 

[71]  -^{TYPE-2  )  /73 

[72]  ->7 4  , CREB*-F3xRCNORM 

[73]  CREB+-  ( F3xRENORMxFCNORM )  tXEMAX 

[74]  ->•  0 

[75]  'WRONG  DATA * 

[76]  +0 

[77]  ’  27/2?  COEFFICIENT  MATRIX  IS  SINGULAR  ’ 

[78]  +  0 
V 


■ 
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